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ABSTRACT 

We report the results of an extensive numerical study of the small-scale turbulent dynamo. The primary focus 
is on the case of large magnetic Prandtl numbers Pr m , which is relevant for hot low-density astrophysical plasmas. 
A Pr m parameter scan is given for the model case of viscosity-dominated (low-Reynolds-number) turbulence. We 
concentrate on three topics: magnetic-energy spectra and saturation levels, the structure of the magnetic-field 
lines, and intermittency of the field-strength distribution. The main results are as follows: (1) the folded structure 
of the field (direction reversals at the resistive scale, field lines curved at the scale of the flow) persists from the 
kinematic to the nonlinear regime; (2) the field distribution is self-similar and appears to be lognormal during the 
kinematic regime and exponential in the saturated state; and (3) the bulk of the magnetic energy is at the resistive 
scale in the kinematic regime and remains there after saturation, although the magnetic -energy spectrum becomes 
much shallower. We propose an analytical model of saturation based on the idea of partial two-dimensionalization 
of the velocity gradients with respect to the local direction of the magnetic folds. The model-predicted saturated 
spectra are in excellent agreement with numerical results. Comparisons with large-Re, moderate-Pr m runs are 
carried out to confirm the relevance of these results and to test heuristic scenarios of dynamo saturation. New 
features at large Re are elongation of the folds in the nonlinear regime from the viscous scale to the box scale and 
the presence of an intermediate nonlinear stage of slower-than-exponential magnetic-energy growth accompanied 
by an increase of the resistive scale and partial suppression of the kinetic-energy spectrum in the inertial range. 
Numerical results for the saturated state do not support scale-by-scale equipartition between magnetic and kinetic 
energies, with a definite excess of magnetic energy at small scales. A physical picture of the saturated state is 
proposed. 

Subject headings: magnetic fields — MHD — plasmas — turbulence — methods: numerical 



1. INTRODUCTION 

1.1. Large- and Small-Scale Magnetic Fields in Astrophysics 

Magnetic fields are detected everywhere in the universe: 
stars, accretion disks, galaxies, and galaxy clusters all carry dy- 
namically important magnetic fields (some of the relevant re- 
views and observations are Kronberg 1994; Beck et al. 1996; 
Minter & Spangler 1996; Zweibel & Heiles 1997; Vallee 1997, 
1998; Balbus & Hawley 1998; Kulsrud 1999; Weiss & To- 
bias 2000; Title 2000; Beck 2000, 2001; Ferriere 2001; Kro- 
nberg et al. 2001; Kronberg 2002; Carilli & Taylor 2002; Han 
& Wielebinski 2002; Widrow 2002; Tobias 2002; Ossendrijver 
2003; Han et al. 2004; Brandenburg & Subramanian 2004). 
There are two kinds of magnetic fields observed. First, there 
are large-scale fields, i.e., fields spatially coherent at scales 
comparable to the size of the astrophysical object that they in- 
habit. Two examples of such fields are the cyclic dipolar field 
of the Sun and the spiral fields of galaxies. Second, there are 
small-scale fields: e.g., the fluctuating fields in the solar pho- 
tosphere, galaxies, and clusters. They are associated with the 



turbulent motions of the constituent plasmas and exist at scales 
below those at which the turbulence is forced. Both types of 
fields are usually strong enough to be dynamically important. 
The challenge is to construct a theory of their origin, evolution, 
and structure consistent with observations and to understand the 
role these fields play in the dynamics of astrophysical objects. 

A physically plausible scenario for the origin and mainte- 
nance of these fields is a turbulent dynamo, which would am- 
plify a weak seed field and culminate in a nonlinear saturated 
state. Just as magnetic fields can be classified into large- and 
small-scale varieties, there are also two kinds of dynamo re- 
sponsible for these fields. First, three-dimensional velocity 
fields sufficiently random in space and/or time will amplify 
small-scale magnetic fluctuations via random stretching of the 
field lines (Batchelor 1950; Zel'dovich et al. 1984; Childress 
& Gilbert 1995). Since it is a small-scale process, we believe 
that this small-scale dynamo can be studied in a homogeneous 
and isotropic setting. In contrast, the generation and structure 
of the large-scale magnetic fields cannot be understood without 
going beyond the homogeneous picture and taking into account 
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large-scale object-specific features: boundary conditions, ro- 
tation, helicity, mean velocity shear, density stratification etc. 
These effects can combine with turbulence to give rise to the 
second kind of dynamo: the large-scale dynamo. This is usu- 
ally handled in the mean-field framework (e.g., Moffatt 1978; 
Parker 1979; Ruzmaikin et al. 1988; Brandenburg & Subrama- 
nian 2004). Mean-field theories tend to assume that all small- 
scale magnetic fluctuations result from the shredding of the 
mean field by the turbulence. Such induced small-scale fields 
do exist, but they vanish if the mean field vanishes and are dis- 
tinct from the dynamo-generated small-scale fields. Neglecting 
the latter is, in fact, not a priori justified because the small-scale 
dynamo is usually much faster than the large-scale one: it am- 
plifies magnetic energy at the rate of turbulent stretching and 
produces dynamically significant fields before the large-scale 
fields can grow appreciably (e.g., Kulsrud & Anderson 1992; 
Kulsrud 1999). Procedures that have been devised for incor- 
porating the effect of the small-scale fields into the mean-field 
theory (Vainshtein & Kichatinov 1983; Radler et al. 2003, and 
references therein) usually do not evolve the small-scale com- 
ponent self-consistently, requiring instead certain statistical in- 
formation about the small-scale fields to be supplied at the mod- 
elling stage. Thus, solving the large-scale problem is predicated 
on making correct assumptions about the small scales. 

Whether the magnetic fields are dynamo-generated or stem 
from primordial and/or external mechanisms (Kulsrud et al. 
1997a; Kronberg et al. 2001), the key question is how they are 
maintained in the course of their interaction with the ambient 
turbulence. Unless the large-scale field is extremely strong, we 
must always ask why it is not quickly churned up by the small- 
scale turbulence. Any large-scale field configuration must be 
consistent with the presence and continued regeneration of con- 
siderable amounts of small-scale magnetic energy. In galaxies 
and clusters, both large- and small-scale fields are in the micro- 
gauss range, which corresponds to magnetic-energy densities 
in approximate equipartition with the kinetic-energy density of 
the turbulent motions in these objects. This suggests that we are 
observing a self-consistent nonlinear state of MHD turbulence. 

While it is surely the interplay between motions and mag- 
netic fields at disparate scales that ultimately determines the 
nature of the observed states, the above considerations motivate 
us to focus on the homogeneous isotropic small-scale MHD tur- 
bulence. From the point of view of theoretical physics, this is 
an attempt to understand the universal features before tackling 
the object-specific ones. From the point of view of numerical 
experimentation, this approach is unavoidable because simulat- 
ing the full picture requires resolving multiple scale separations 
that are well beyond the capabilities of current computers and 
are likely to remain so for many years to come. 

1.2. The Scales in the Problem 

The issue of multiple scale ranges is important, so let us ex- 
amine it in more detail. The hierarchy of scales in an astrophys- 
ical object can be outlined as follows. The largest spatial scale 
is the system size L. The energy that feeds the turbulence is 
injected at the outer, or forcing, scale Iq. For objects with non- 
trivial geometry such as stars or galactic disks, Iq <C L. This is 
the first scale separation in the problem. It is used by mean-field 
theories to split all fields into mean and fluctuating components 
with averages being done over scales of order and below Iq. 

The energy cascades from Iq down to the viscous-dissipation 
scale l v . In purely hydrodynamic systems, the latter is called 
the Kolmogorov inner scale. Kolmogorov's 1941 dimensional 



theory of turbulence (e.g., Frisch 1995) gives /„ ~ Re~ 3//4 Zo- The 
hydrodynamic Reynolds number Re ~ uqIq /v is usually fairly 
large for astrophysical systems, so l u -C Iq (uq is the typical ve- 
locity at the outer scale, v is the fluid viscosity). This is the 
second scale separation in the problem. It is at scales between 
the outer and the inner scale — the inertial range — that the 
universal physics of turbulence is contained. 

If magnetic fields are present, they bring in their own dis- 
sipation scale L associated with the magnetic diffusivity ij of 
the plasma (Spitzer 1962). The ratio Pr m = v/rj is called the 
magnetic Prandtl number. Using the Spitzer values for the vis- 
cosity and magnetic diffusivity of a fully ionized plasma, one 
finds Pr m ~ 10~ 5 r 4 /n, where T is temperature in K and n is 
the particle concentration in cm" 3 . In the case of a partially 
ionized gas with neutral-dominated viscosity (e.g., warm inter- 
stellar medium), we can estimate v ~ v^/an, where v t h is the 
thermal speed and a the atomic cross-section. The formula for 
Pr m is then Pr m ~ 10 7 T 2 /n. In hot low-density plasmas such as 
(warm) interstellar and intracluster media, as well as in some 
accretion disks, Pr m 3> 1 (see, e.g., Brandenburg & Subrama- 
nian 2004). Systems with Pr m 3> 1 will be the primary focus 
of this paper. If a weak magnetic field is amplified by Kol- 
mogorov turbulence, the magnetic-energy exponentiation time 
is the same as the turnover time of the viscous-scale eddies, be- 
cause these eddies are the fastest ones. Balancing the viscous- 
eddy turnover time with the magnetic-field diffusion time gives 
an estimate for the resistive scale in the kinematic (weak-field) 
regime: l n ~ Pr" 1 / 2 /^. As Pr m 3> 1, l v <C l v . This is the third 
scale separation in the problem. The scale range in between 
contains degrees of freedom that are accessible to magnetic 
fluctuations but not to velocities (Fig. 1). It is hardly surprising 
that these degrees of freedom are quickly occupied: the small- 
scale kinematic dynamo spreads magnetic energy over the sub- 
viscous range and piles it up at the resistive scale (e.g., Kulsrud 
& Anderson 1992). The key question then is what happens 
when the small-scale fields become sufficiently strong for the 
Lorentz force to react back on the flow. 
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FIG. 1. — Sketch of scale ranges and energy spectra in a large-Pr m medium. 

We have identified three scale ranges, L ^> Iq ^> l v ^> l n . To 
give a concrete example, let us give order-of-magnitude esti- 
mates for these scales in our Galaxy. The disk diameter is 
L ~ 10 4 pc; the supernova scale at which the turbulence is 
stirred is Iq ~ 10 2 pc; the Reynolds number is Re ~ 10 5 , so the 
viscous scale is l v ~ 10~ 2 pc; the Prandtl number is Pr m ~ 10 14 , 
so the resistive scale is 10 4 km, a tiny distance by Galactic stan- 
dards (note §2.6). The associated timescales are the period 
of Galactic rotation T ~ 10 s yr, the outer-eddy turnover time 
To ~ 10 7 yr, and the viscous-eddy turnover time t v ~ 10 5 yr. 
The mean-field dynamo theory gives a large-scale Galactic field 
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exponentiating in a rotation time T, while the timescale for the 
growth of the small-scale magnetic energy is t„ (Kulsrud & An- 
derson 1992). This illustrates our earlier point that the small- 
scale dynamo tends to be faster than the large-scale one. 

1.3. Simulating the Problem 

Simultaneously resolving all of these scales is not achievable 
in numerical simulations. Even if the large-scale object-specific 
features are forsaken and only isotropic homogeneous MHD 
turbulence with large Re and Pr m is considered in a periodic 
box, it is still not possible to simulate both the hydrodynamic 
inertial and the magnetic subviscous ranges in the same box. 

In the relevant numerical work so far, a popular course of 
action has been to run simulations with Pr m = 1 (Meneguzzi et 
al. 1981; Kida et al. 1991; Kleva & Drake 1995; Miller et al. 
1996; Cho & Vishniac 2000; Mininni et al. 2003; Haugen et 
al. 2003, our runs al, a2, and A in § 5). This choice is usually 
motivated by the view of MHD turbulence as a cascade of in- 
teracting Alfven-wave packets (Kraichnan 1965; Goldreich & 
Sridhar 1995). In an Alfven wave, velocity and magnetic fluc- 
tuations are equal, dissipation occurs at either viscous or resis- 
tive scale, whichever is larger, and whatever happens below that 
scale is usually not expected to affect the inertial-range physics. 

In fact, this theory only appears to work for anisotropic 
MHD turbulence with an externally imposed strong mean field 
(Maron & Goldreich 2001 ; Cho et al. 2002; Muller et al. 2003). 
In (forced) simulations with zero mean field, one invariably sees 
that velocity and magnetic field are not symmetric, their statis- 
tics are not the same, they do not dissipate at the same rate or 
at the same scale, and there is no scale-by-scale equipartition 
of kinetic and magnetic energies, as there should have been if 
the turbulence had been purely Alfvenic (Maron et al. 2004, and 
see §5.1 .2). The picture that does emerge is rather that of a very 
significant excess of magnetic energy at small scales and the 
magnetic-energy spectra resembling the resistively-dominated 
spectra of the small-scale large-Pr m kinematic dynamo. 

A natural step beyond Pr m = 1 is to look at Prandtl numbers 
above, but of the order of, unity (Meneguzzi et al. 1981; Chou 
2001; Brandenburg 2001; Maron & Blackman 2002; Maron 
et al. 2004; Archontis et al. 2003a,b; Haugen et al. 2004, our 
run B in § 5) more or less in the hope that essential features of 
the Pr m ^ 1 regime would be captured. This way, the inertial 
range might be resolvable, but the subviscous one (whose width 
is <~ PrJ/ 2 ) is largely sacrificed. This approach is, to some ex- 
tent, justified (see § 5). However, as the asymmetry between 
magnetic and velocity field becomes more pronounced with in- 
creasing Pr m , one is faced with the realization that the Pr m > 1 
case represents a nonasymptotic mixed state about which it is 
very hard to make any clear statements. We have pointed out 
in an earlier paper (Schekochihin et al. 2002d, see also § 5.2.3) 
that, in order for the effects of the scale separation between l v 
and l v on the nonlinear physics to be fully present, one must 



have Pr m 3> Re 1 / 2 >• 1. Astrophysical plasmas have no prob- 
lem satisfying this condition, but numerically simulating this 
regime is not, as yet, possible. 

We believe that the key feature of isotropic MHD turbulence 
is the scale separation that arises between velocities and mag- 
netic fields. In an atrophied nonasymptotic form, this scale sep- 
aration is discernible even at Pr m = 1 . Therefore, in most of this 
paper (with the exception of § 5), we follow Kinney et al. (2000) 
in choosing to resolve the subviscous range while sacrificing the 
inertial one. Namely, we set up our simulation in such a way 
that the forcing and viscous scales are comparable, while most 
of the available resolution is spent on the subviscous range. 5 At 
resolutions of 256 3 , we can study Prandtl numbers up to 2500, 
which just barely allows us to make some statements about the 
asymptotic Pr m 1 limit. Furthermore, this approach allows 
us to study Pr m dependence of our model, which is the only pa- 
rameter scan that can at present be afforded. Some preliminary 
results obtained in the same setting (at lower resolutions) were 
published in Schekochihin et al. (2002c). 

1.4. Plan of the Paper 

Our numerical set up, equations, and the limitations of the 
model are discussed in more detail in § 2. Section 3 is de- 
voted to a detailed study of the kinematic regime, which will be 
important for understanding further developments. Section 4 
treats the nonlinear saturated state. Summaries are provided at 
the end of each of these two sections to help the reader identify 
the key points. In § 5, our understanding of the physics of the 
nonlinear dynamo with large Reynolds numbers is outlined and 
some tentative comparisons are made with the results of simu- 
lations with large Re and Pr m of order unity. Section 6 discusses 
the scope of applicability of the large-Pr m -dynamo results and 
discusses directions of future work. An itemized summary of 
the main results and conclusions of the paper is given in § 7. 

2. THE MODEL 

2. 1 . The Equations and the Code 

We consider the equations of incompressible MHD: 

d_ 

dt 
d_ 

dt ' 

where d /dt = d, + u • V is the advective time derivative, u(f , x) is 
the velocity field, B(f ,x) is the magnetic field, and f(f ,x) is the 
body force. The pressure gradient Wp (which includes mag- 
netic pressure) is determined by the incompressibility condi- 
tion V • u = 0. We have normalized p and B by p and (Airp) 1 / 2 , 
respectively, where p = const is plasma density. 

Equations (1-2) are solved in a triply periodic box by the 
pseudospectral method (see code description in Maron & Gol- 
dreich 2001; Maron et al. 2004). The body force f is random, 



: ^Au-V/? + B- VB + f, 
B=B V11 + 77AB, 



(1) 



(2) 



5 Note that another approach in which low Reynolds numbers are often considered is to study magnetic-field growth and saturation in, and nonlinear modifications 
to, specific deterministic time-independent or time-periodic velocity fields with chaotic trajectories that are known to be dynamos in the kinematic limit (Childress & 
Gilbert 1995). The best-documented velocity fields of this sort are the ABC flows. The strategy then is to set up a forcing function that, in the absence of the magnetic 
field, would reproduce a given flow as a solution to the Navier-Stokes equation, and to keep Re below the stability threshold of this solution — or to use some 
simplified model of the velocity equation that accomplishes the same, — and to model the fast-dynamo regime by maximizing Pr m (Galanti et al. 1992; Cattaneo et al. 
1996; Maksymczuk & Gilbert 1999). This approach differs from ours in that the velocity field is not random in time, nor is it homogeneous in space, so such features 
as stagnation points of the flow, transport barriers, etc., may acquire a degree of importance that cannot be a feature of real turbulence. The redeeming advantage here 
is in being able to deal with a velocity field that, unlike in the case of real turbulence, does not constitute an unsolved problem by itself, so one has more control over 
the numerical experiment and some of the nonlinear effects may be more straightforward to detect. In some of these studies, Re is allowed to exceed the hydrodynamic 
stability threshold (Zienicke et al. 1998; Brummell et al. 2001), in which case the small-scale results are probably universal because of the usual physical expectation 
that the small-scale properties of turbulence are insensitive to the particular choice of large-scale forcing. Of course, in order for such a universality argument to work, 
Re must be sufficiently large, and resolving large Prandtl numbers is again problematic. 



4 



SCHEKOCHIHIN ET AL. 



nonhelical, applied at the largest scales in the box (ko/2-K = 
1,2), and white in time (i.e., statistically independent at each 
time step): 

(f(t,x)f(t',x')) = <5(f-f')e' 7 (x-x'). (3) 
For a white-noise forcing, the average injected power is fixed: 

^(" 2 )=-^(|Vu| 2 )-(BB:Vu) + 6, (4) 

where e = (u • f) = (1 /2)e"(0) = const. The code units are based 
on box size 1 and injected power e = 1 . 

All runs are listed in Table 1 . Runs S0-S6 are viscosity domi- 
nated (§ 2.3) with S4 being the main time-history run. Runs al, 
a2, A, and B have larger Re but smaller Pr m (see § 5). 

2.2. Averaging 

Mathematically speaking, all fields are random because the 
forcing is random. The system is isotropic and homogeneous in 
space, so all vector fields have zero mean and all points in the 
box are equally nonspecial. In theory, the angle brackets mean 
ensemble averages with respect to the forcing realizations and 
to the initial conditions. In practical measurements, one can 
(based on the usual ergodic assumption) use volume and time 
averaging instead. In what follows, wherever time histories are 
plotted, only volume averaging is done, while for quantities av- 
eraged over the kinematic or saturated stages of our runs, both 
volume and time averaging are performed. Error bars show the 
mean square deviations of volume averages at particular times 
from the time-averaged values. 

Averages in the kinematic regime are over 20 code time units 
(51 snapshots with time separation Af = 0.4) for runs S0-S4, 
over 6 time units (31 snapshots, At = 0.2) for runs A and B, 
over 35 time units (71 snapshots, Af = 0.5) for run al, and 
over 15 time units (31 snapshots, Af = 0.5) for run a2; aver- 
ages in the saturated state are over 20 time units (5 1 snapshots, 
Af = 0.4) for runs S1-S6, over 10 time units (51 snapshots, 
Af = 0.2) for runs A and B, and over 30 time units (61 snap- 
shots, Af = 0.5) for runs al and a2. We have checked that in- 
creasing the averaging intervals does not alter the results. One 
code time unit roughly corresponds to the box crossing time 
(roughly the turnover time of the outer eddies). 

2.3. The Viscosity-Dominated Limit 

As we can only adequately resolve either the inertial or the 
subviscous range, but not both, we choose to concentrate on the 
latter. In the runs reported in § 3 and § 4 (but not in § 5), the vis- 
cosity is fixed at v = 5 x 10" 2 , which effectively leads to lo~l v . 
The velocity field in these runs is not, strictly speaking, turbu- 
lent. It is still random in time because the external forcing is 
random, but it is smooth in space as a result of very strong vis- 
cous dissipation. Advection by a random but spatially smooth 
flow is known as the Batchelor regime (Batchelor 1959). From 
the point of view of realistic turbulence, such a flow is probably 
a fairly good model of the viscous-scale eddies acting on the 
magnetic fluctuations in the subviscous range. In this interpre- 
tation, the external forcing f models the energy supply from the 
Kolmogorov cascade to the motions at the viscous scale. 

We emphasize that, while f is white in time, the resulting ve- 
locity is not. In order to see this, let us consider eq. (1) in the 
more drastic limit of Re <C 1, so the hydrodynamic nonlinear- 
ity can be neglected, and let the Lorentz forces be negligible as 
well. Then the solution of eq. (1) in k space is 

u(f,k) = e- l/ * 2 'u(f = 0,k)+ / dTe-" k2r f(t-T,k). (5) 
Jo 



If f is Gaussian, so is u, but u has a finite correlation time: 

( M '(f,k)«>(f',k')) = (27r) 3 (5(k+k')^ e -^ 2 "-''le^(k), (6) 

where e y (k) = J d 3 ye~ lky e'j(y) is the Fourier transform of the 
forcing correlation tensor [eq. (3)], and we tookf ,f' 3> {vk,Q 2 )~ x , 
Icq being the forcing wavenumber. The velocity correlation time 
is then r c ~ {vk 2 )~ l . Equations (5-6) are strictly valid only if 
vk 2 ^> kuk for all k. In our simulations, we chose v so that 
vko 2 ~ koUk (•<=> k v <~ ko), i.e., just large enough for the inertial 
range to collapse. The obvious estimate for the velocity corre- 
lation time in this regime is t c ~ ikau^Y 1 . Given the limitation 
of not having an inertial range, we believe this to be a sensible 
way of modelling turbulent motions. 

Explicit contact can be made between the viscosity- 
dominated model and the analytical results from the kinematic- 
dynamo theory. One of the very few time-dependent random 
flows for which kinematic dynamo is exactly solvable is a Gaus- 
sian white velocity field t)(f ,x): 

(t3\t,\)X) j (t',x')) = S(t-t')K ij (x-x'). (7) 

This model was introduced by Kazantsev (1968) and, in the 
context of passive-scalar advection, independently by Kraich- 
nan (1968). For the Pr m 1 problems, many quantities of in- 
terest then turn out to depend only on the first few coefficients 
in the Taylor expansion of the velocity correlation tensor: 

« y (y)= ^-\K 2 {y 2 8^-y y i) 

i^mv'(v'^ T.v'rh... (8) 

Increasing v in eqs. (5-6) above the value at which k v ~ ko is 
equivalent to decreasing the correlation time of the flow. Thus, 
the velocity field in the viscosity-dominated model reduces to 
the Kazantsev velocity in the limit v — > oo, provided that the 
force is rescaled to keep the integral of the time-correlation 
function constant: e !J (k) = v 2 k A n l '()^) (i.e., f = -i>Ao). Since 
we owe much of our understanding of the kinematic dynamo 
to the Kazantsev model (e.g., Kazantsev 1968; Artamonova 
& Sokoloff 1986; Kulsrud & Anderson 1992; Gruzinov et 
al. 1996; Chertkov et al. 1999; Schekochihin et al. 2002b,d, 
2004b), it is convenient to be able to consider the nonlinear 
problem as a well-defined extension of it. 



2.4. Incompressibility 

We use incompressible MHD even though most astrophys- 
ical flows are not incompressible at large scales. For exam- 
ple, the supernova forcing of the galactic turbulence produces 
sonic velocities at the outer scale. However, motions at smaller 
scales are subsonic and mostly vortical and incompressible 
(e.g., Boldyrev et al. 2002; Porter et al. 2002; Balsara et al. 
2004). As we are interested in the small-scale dynamo, which is 
driven by the smallest eddies, the incompressibility assumption 
is, therefore, reasonable. This view is further supported by the 
satisfactory outcome of a number of comparisons (carried out 
by N. E. Haugen 2003, private communication) between our 
simulations (§ 5) and the simulations of Haugen et al. (2003, 
2004), who used a compressible grid code. 
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2.5. Helicity 

Our simulations are nonhelical in the sense that the net ki- 
netic helicity is zero, (u • (V x u)) = (helicity fluctuations are 
allowed). The net helicity of the flow can be controlled via the 
helical component of the external forcing (Maron & Blackman 
2002). Helicity is usually considered important because it is 
present in some of the astrophysical objects of interest. How- 
ever, net helicity is a large-scale feature associated with overall 
rotation. Its effects should only be felt at timescales compara- 
ble to the rotation timescale, which is usually much longer than 
the turnover time of the turbulent eddies. Furthermore, helicity 
undergoes an inverse cascade from small scales to the largest 
scale in the box (Frisch et al. 1975; Pouquet et al. 1976). While 
net helicity is crucial for the generation of the box-scale mag- 
netic field (essentially the mean field), the above considerations 
suggest that it probably does not play a significant role in the 
small-scale dynamo. In the viscosity-dominated simulations re- 
ported in § 3 and § 4, the box scale is the viscous scale and the 
external forcing models energy influx from larger scales, so al- 
lowing helical forcing is not a physically sensible choice. In the 
large-Re runs of § 5, the effect of helicity would be interesting 
to study, but it would introduce another scale separation into 
the problem (between the box size and the forcing scale) and, 
therefore, require dramatic increases in numerical resolution. 

2.6. Plasma Dissipation Processes 

Our choice of Laplacian diffusion operators to describe vis- 
cous and resistive cutoffs is only a very crude model of the ac- 
tual dissipation processes in astrophysical plasmas. 

In a fully-ionized plasma, ions become magnetized already 
at relatively low magnetic-field strengths and the hydrodynamic 
viscosity must be replaced by a more general tensor viscosity 
(Braginskii 1965, see also Montgomery 1992), which is locally 
anisotropic and very inefficient at damping velocity gradients 
perpendicular to the magnetic field: the perpendicular compo- 
nents of the viscosity tensor vanish in the limit of zero Larmor 
radius (the magnetic diffusivity also becomes anisotropic, but 
the difference between its perpendicular and parallel compo- 
nents is only a factor of <~ 2). This can alter the physics at 
subviscous scales (Malyshkin & Kulsrud 2002; Balbus 2004). 

In partially ionized plasmas (e.g., in the interstellar medium), 
viscous dissipation is controlled by neutral species and can be 
assumed isotropic. However, in the presence of neutrals, am- 
bipolar damping might enhance magnetic-field diffusion (e.g., 
Zweibel 1988; Kulsrud & Anderson 1992; Proctor & Zweibel 
1992; Brandenburg & Zweibel 1994; Brandenburg & Subrama- 
nian 2000; Zweibel 2002; Kim & Diamond 2002). 

Finally, since the mean free path in astrophysical plasmas is 
often much larger than the scale associated with the Spitzer re- 
sistivity, the magnetic-field diffusion could be superseded by 
kinetic damping effects not contained in the MHD description 
(Foote & Kulsrud 1979; Kulsrud et al. 1997b). 

These and other effects outside the MHD paradigm can be 
important. However, in this paper, we have opted to study the 
problem in the simplest available formulation: that provided by 
eqs. (1-2). We believe that this minimal model already contains 
the essentials of the small-scale dynamo physics. 

3. THE KINEMATIC DYNAMO 

The kinematic dynamo is a natural starting point for a nu- 
merical experiment. What we learn about magnetic fields in the 
kinematic limit will motivate our investigation of the nonlinear 



regime. This is a rather detailed study, so a busy reader anxious 
to get to the nonlinear matters may skip to the short summary 
in § 3.5 and on to § 4. 




Straight field Folded field 



FIG. 2. — Stretching and folding of field lines by turbulent eddies. 

The small-scale turbulent dynamo is caused by the random 
stretching of the (nearly) frozen-in field lines by the ambient 
random flow. Stretching leads to exponential growth of the field 
strength (Batchelor 1950). The growth rate 7 ~ Vu is then nat- 
urally of the order of the inverse eddy-turnover time. In the 
case of Kolmogorov turbulence, it is the turnover time t v of 
the smallest (viscous-scale) eddies, because they are the fastest 
(see §5.1.1). For the magnetic fields at subviscous scales, the 
action of these eddies is roughly equivalent to the application 
of a random linear shear. Without diffusion (r\ = 0), stretch- 
ing is not opposed by any dissipation mechanism. However, 
stretching is necessarily accompanied by the refinement of the 
field scale (Fig. 2), which proceeds exponentially fast in time 
and also at the eddy-turnover rate. Thus, the field scale even- 
tually becomes comparable to the diffusion (resistive) scale, 
bringing an end to the diffusion-free evolution. The charac- 
teristic time during which diffusion-free considerations apply 

is t - 7 " 1 ln(^Ar,)^7" 1 ln(Pri 1 /2 ), assuming the initial field 
varies at the velocity scales. 
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FIG. 3. — Exponential growth and saturation of the magnetic energy. Full 
time history is only shown for run S4. See Table 1 for the index of runs. 

Once diffusion enters the picture, it will (partially) offset the 
effect of random stretching (see § 3.3). If magnetic energy con- 
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FIG. 4. — (a) Evolution of the magnetic-energy spectrum for run S4. The spectra for the diffusion-free regime (0 < t < 2) are given at time intervals of Af = 0.4. 
The subsequent evolution (2 < t < 40) is represented by spectra at time intervals At = 2. (b) Magnetic-energy spectra (normalized by the total magnetic energy and 
averaged) during the kinematic stage. Note that for run SO dynamo is resistively suppressed. 



tinues to grow exponentially in this regime with a growth rate 
that does not vanish as r\ — > +0, the dynamo is called fast (Vain- 
shtein & Zel'dovich 1972). Whether any given flow is a fast 
dynamo is a problem that very rarely has an analytical solu- 
tion (Childress & Gilbert 1995). In practice, since the advent 
of numerical simulations, turbulent flows have more often than 
not been found to support fast dynamo action (see references 
in § 1.3 and discussion in § 6.1). In what follows, we shall see 
that magnetic fluctuations that result from this dynamo have a 
distinctive structure (see § 3.2 and § 3.3). 

In all our runs, the initial magnetic field had a spectrum con- 
centrated at velocity scales so that both diffusion-free and diffu- 
sive regimes of the kinematic dynamo are realized. As a result 
of resolution constraints, the diffusion-free stage is never longer 
than 2 time units. The diffusive stage is quite long, so it is pos- 
sible to extract statistical information by averaging over 20 time 
units. The magnetic energy grows exponentially from the initial 
value of ~ 10" 10 to saturated values of order unity (Fig. 3). 

In this Section, we study the properties of the growing field. 
We concentrate on three main questions. First, what is its en- 
ergy spectrum? This is considered in § 3.1. Second, what does 
the field "look" like? The field structure is studied in § 3.2. 
Third, how does it fill the volume? This question is treated in 
§ 3.4 in terms of the field-strength statistics. § 3.3 is a qualita- 
tive interlude discussing the physical reasons for the existence 
of the small-scale dynamo. 



3.1. The Magnetic-Energy Spectrum and Growth Rate 

The evolution of the angle-integrated magnetic-energy spec- 
trum, M(k) = (1 /2) Jdn k k z (\B(k)\ 2 }, in the large-Pr m regime is 
the best understood part of the kinematic dynamo physics. For 
the Kazantsev velocity model, M(k) can be shown to satisfy a 
closed integrodifferential equation, which, for k ^> k u , reduces 
to a simple Fokker-Planck equation (Kazantsev 1968; Vain- 
shtein 1972, 1980, 1982; Kulsrud & Anderson 1992; Gruzinov 



et al. 1996; Schekochihin et al. 2002a) 



d,M= - — k 2 — M-4-kM + 2^M-2nk 2 M, (9) 
5 ok \ ok ) 

where 7 = (5 /4)k2 ~ Vu is comparable to the inverse turnover 
time of the viscous eddies. Integrating eq. (9), gives 

d l (B 2 )=2j(B 2 )-2 V k 2 m ,(B 2 ), (10) 

so 7 is the growth rate of the rms magnetic field in the ab- 
sence of diffusion (the stretching rate). Without diffusion, 
eq. (9) describes a spreading lognormal profile with a peak 
moving exponentially fast toward ever larger wavenumbers, 
& P eak oc exp[(3/5)7?], the amplitude of each mode growing ex- 
ponentially at the rate (3/4)7, and the dispersion of In k increas- 
ing linearly in time as (4/5)7?. 

Once the resistive scale is reached and diffusion becomes 
important, further increase of k stops. The spectrum is a solu- 
tion of an eigenvalue problem for eq. (9) with M(k) oc exp(A7? )• 
Here A < 2 because diffusion dissipates some of the dynamo- 
generated energy: eq. (10) is an expression of this competition 
between stretching and diffusion. If A tends to a positive r\- 
independent constant in the limit r\ — > +0, stretching wins and 
the dynamo is fast. In order to solve this eigenvalue problem, 
we must specify the boundary condition at small k. In princi- 
ple, this boundary condition can be determined only by solving 
the original integrodifferential equation of which eq. (9) is an 
asymptotic form for » This integrodifferential equation 
contains the velocity statistics at all scales and cannot be solved 
in a general and model-independent way. However, the eigen- 
values are not very sensitive to the boundary condition: a nat- 
ural choice is to impose a zero-flux constraint at some infrared 
cutoff k* (Schekochihin et al. 2002d), but a number of other 
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reasonable options give the same solution 6 

M(k) ~ (const)e A ^fe 3/2 Xo(fe/^), 
3 
4 



A~ -- 



(11) 
(12) 



short-scale correlations across. The simplest approach is, in 
fact, to measure average characteristic scales at which the field 
varies in the directions perpendicular and parallel to itself. In 
two dimensions, this was done by Kinney et al. (2000). In three 
dimensions, the characteristic scales are studied in §3.2.1. A 
more detailed description of the field's "statistical geometry" is 
provided by the distribution of the field-line curvature (§ 3.2.4) 
and by its relation to the field strength (§ 3.2.3). The advantage 
of this approach is that curvature is a local quantity, so we only 
have to look at one-point statistics. It is also directly involved in 
the Lorentz tension force [eq. (27)] thus providing information 
about the onset of the nonlinear back-reaction (§ 4). An analyt- 
ical theory based on the Kazantsev model was developed in an 
earlier paper (Schekochihin et al. 2002b). The numerical study 
below is motivated by the quantitative predictions made there. 



V 



5[ln(V**)] 2 ' 

where k v = (7/10?7) 1//2 and is the Macdonald function. In the 
limit rj — > +0, we get A — > 3/4, a fast dynamo. Note that (3/4)7 
is the rate at which individual £-space modes grow as a result 
of random stretching, so the effect of diffusion is simply to stop 
the magnetic energy from spreading to ever larger k (Kulsrud & 
Anderson 1992; Schekochihin et al. 2002a). 

The numerical results broadly confirm this theoretical pic- 
ture. In Figure 3, we show the evolution of the magnetic en- 
ergy for several values of Pr m . The growth is exponential in 
time and, for Pr m = 500, the growth rate after resistive scales 
are reached (f > 2) is reduced compared to the diffusion-free 
regime by a factor that is actually quite close to 3/8 predicted 
by the Kazantsev model. 7 The growth rate at smaller values of 
Pr m is smaller, as a result of nonnegligible reduction by diffu- 
sion as suggested by the logarithmic correction in eq. (12). 

The Prandtl number must exceed a certain critical value Pr m c 
in order for the dynamo to be possible. Indeed, in eq. (12), for 
k* ~ k v , we have k v /k* ~ PrJ/ 2 and setting A = gives a rough 
estimate of the threshold below which the dynamo is resistively 
suppressed: Pr m c ~ 26. In our runs, we have found growth 
at Pr m = 50 and decay at Pr m = 25, so Pr mc e (25,50). Note 
that both the Kazantsev velocity with smooth correlation ten- 
sor (8) and the velocity field in our simulations are single-scale 
flows with Re ~ 1, so Pr m c simply corresponds to the criti- 
cal magnetic Reynolds number Rm c = RePr m c . Since for our 
viscosity-dominated runs Re ~ 2, we have Rm c 6 (50, 100). 

The diffusion-free spreading of the magnetic energy toward 
resistive scales and the subsequent self-similar growth of the 
spectrum with a peak at the resistive scale [eq. (11)] are man- 
ifest in Figure 4a. The k^l 2 scaling appears to be correct 
(Fig. 4b), although resolving a broader scaling range is nec- 
essary for definite corroboration. The occasional strong disrup- 
tions of the self-similar growth (Fig. 4a) are discussed at the 
end of §3.4.3. 

3.2. The Folded Structure 

The small-scale magnetic fields produced by the dynamo are 
not, in fact, completely devoid of large-scale coherence. The 
smallness of their characteristic scale is due primarily to the 
rapid (in space) field reversals transverse to the local field di- 
rection. Along themselves, the fields vary at scales comparable 
to the size of the eddies. This folded structure is schemati- 
cally illustrated in Figure 2 and is evident in field visualizations 
(Fig. 5). 

There are many ways to diagnose the folded structure. Ott 
and coworkers studied the field reversals in map dynamos in 
terms of magnetic-flux cancellation: see review by Ott (1998) 
and numerical results by Cattaneo et al. (1995). Chertkov et al. 
(1999) considered two-point correlation functions of the mag- 
netic field and found large-scale correlations along the field and 

6 These include, e.g, requiring that the magnetic-field second-order correlation function decay exponentially at large distances (Artamonova & Sokoloff 1986; Sub- 
ramanian 1997) or setting M(k„) = (Schekochihin et al. 2002a). Kulsrud & Anderson (1992) solved the integrodifferential equation for the spectrum numerically, 
using a Kolmogorov spectrum for velocity, and also obtained A = 3/4. See also Gruzinov et al. (1996) for yet another argument leading to the same result. 

7 We hesitate to claim that this is a corroboration of the specific number that obtains in the Kazantsev model rather than a mere coincidence. The Kazantsev velocity 
is (^-correlated in time, which cannot be a good quantitative description of the velocity field in our simulations (see § 2.3). It is well known that finite-correlation-time 
effects can result in order-one corrections to the dynamo growth rate (Chandran 1997; Kinney et al. 2000; Schekochihin & Kulsrud 2001; Chou 2001; Kleeorin et al. 
2002). While the scalings in the Kazantsev model may be universal (e.g., the power tail of the curvature pdf, § 3.2), the values of growth rates can only a priori be 
considered as qualitative predictions. 



FIG. 5. — Cross section of the magnetic field strength (gray scale) and in- 
plane field direction (red arrows) at t = 20 during the kinematic stage of run S4. 
Lighter regions correspond to stronger fields. Note that what appears to be 
strong-field "clumps" are, in fact, cross sections of folds transverse to the page. 

3.2.1. Characteristic Scales 
We define the characteristic parallel wavenumber of the field 



Ml 



IB • VBI 



(B 4 ) 



1/2 



and its characteristic perpendicular wavenumber 

; o > .11 

^BxJ : 



(B 4 ) 



(13) 



(14) 



where J = V x B. The overall field variation is measured by the 
rms wavenumber 



1 



(B 2 )/2 



dkk 2 M(k) 



1/2 



IVBI 
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(15) 
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FIG. 6. — (a) Evolution of characteristic wavenumbers (defined in § 3.2.1) and of AT rms = (|b • Vb| 2 }'/ 2 for ran S4. (b) Averaged values of the same wavenumbers 
vs. Pr m during the kinematic stage of runs S0-S4. Note that for run SO dynamo is resistively suppressed. The values plotted are listed in Table 1. 

We find (Fig. 7) that all these spectra grow self-similarly and 



While both £bxj an d ^™ grow exponentially during the 
diffusion-free regime and saturate at values that scale as k n ~ 
fc^Pr^ 2 (cf. Brummell et al. 2001), the parallel wavenum- 
ber fen remains approximately constant and comparable to k v 
(Schekochihin et al. 2002b): see Fig. 6. For comparison, we 
also plot (Fig. 6) the rms wavenumber of the velocity field, 



1 



<«2)/2 



dkk 2 E(k) 



1/2 



(|Vu| 2 ) 

<« 2 > 



1/2 



V5 
A 



(16) 



where E(k) is the angle-integrated kinetic-energy spectrum 
and A is the Taylor microscale of the flow defined in the 
standard way (Frisch 1995). In Kolmogorov turbulence, A <~ 
~ Re 1/4 l v , so k x ~ Re^X- 
In order to estimate the degree of misalignment between the 
direction-reversing fields, it is instructive to look at a character- 
istic wavenumber of field variation in the direction orthogonal 
to both B and B x J: 



J = 



|B-J| 2 ) 

(B 4 ) 



-,1/2 



(17) 



Figure 6 and Table 1 show that k^.j <~ k u and does not increase 
with Pr m in the kinematic regime. Thus, the reversing straight 
fields are fairly well aligned and the flux surfaces are sheets. 

3.2.2. Fourth-Order Spectra 

One can probe further into the behavior of the quantities 
B • VB, B x J, and B J by looking at their spectra. All these 
spectra are, in fact, controlled by the spectrum of B 2 , which 
turns out to be flat. Define angle-averaged spectra 



M 4 (k) ~ k°, M BxJ (k) ~ i k 2 M A (k), (22) 

M B .j(k) ~ kl. } M 4 (k), T(k) ~ k 2 M 4 (k), (23) 

with and ^b j defined by eqs. (13) and (17). 

This simple behavior is due to the folded structure. If the 
magnetic-field variation is dominated by direction reversals, 
then, for k v <C k <C k n , the main contribution to the integral 

(|[B 2 ](k)| 2 ) = 

J d 3 k' J d 3 k" (B(k')-B(k-k')B(k")-B(-k-k"))(24) 

is from k! ,lc" ~ k v with the proviso that k', k", and k are all in 
the direction of the reversal (i.e. transverse to the flux sheet). 
Expanding in k, we get to zeroth order, (|[B 2 ](k)| 2 ) ~ (B 4 ) cx 
k°. Since the spectrum is essentially one-dimensional (peaked 
for k transverse to the flux sheet), we get M 4 (k) <~ k°. The rest 
of the fourth-order spectra can be worked out in the same way 
and eqs. (22-23) are readily obtained. 

These results prove useful in § 4.2 and § 5.3.4. 

3.2.3. Magnetic-Field Strength and Curvature 

The field geometry can be studied in a more detailed way in 
terms of the curvature K = b Vb of the field lines (b = B/B). 
K= |K| satisfies (cf. Drummond & Munch 1991) 





1 dn k k 2 (\[B 2 ](k)\ 2 ), 


(18) 


Afnxj(*)=£ 


fdn k k 2 (\[BxJ](K)\ 2 }, 


(19) 


M*.m = \ 


f ^ k fe 2 (|[B-J](k)| 2 ), 


(20) 


T(k)= 1 - 


f ^ k fe 2 (|[B-VB](k)| 2 >. 


(21) 



K= (nn: Vu-2bb: Vu)7r+bb: (VVu) 



(25) 



dt 

where n = K/K, and resistive terms have been dropped. Equa- 
tion (25) is a straightforward consequence of eq. (2). A stan- 

9 1/2 

dard kinematic calculation of the rms curvature, K [ms = (K ) , 
shows that it grows exponentially in the diffusion-free regime 
(Malyshkin 2001; Schekochihin et al. 2002b): see Fig. 6. Since 
curvature is an inverse scale, it is clear that in the diffusive 
regime, K rms must saturate at some Pr m -dependent value. 
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FIG. 7. — (a) Tension spectra for run S4 for 8 < t < 40 at intervals At = 2 and in the saturated state, (b) Fourth-order spectra defined in eqs. (18-21) (normalized 
by (B 4 )/2 and averaged) during the kinematic stage of run S4. The evolution of all these spectra is analogous to that of T(k) shown in (a). 



The growth of K lms does not contradict our previous claim 
that magnetic field lines are straight up to the flow scale. As the 
field is stretched and folded, it is organized in long thin struc- 
tures (folds) where it is only significantly curved in the bends 
(turning points). The flow stretches the straight segments of 
the field, while the curved fields in the bends remain weak (this 
is simply a consequence of flux and volume conservation: see 
Schekochihin et al. 2002b, c). Thus, there is anticorrelation be- 
tween the field-line curvature and the field strength. 

The most straightforward statistical measure of this anticor- 
relation is the correlation coefficient defined as follows 
(K 2 B 2 ) - (K 2 )(B 2 ) 



is quite wide. It is clear, however, that magnetic fields with 
B > firms have curvatures well below 7T rms while the fields with 
curvatures K > K lms are quite weak. 

If we consider eq. (25) together with the evolution equation 
for the magnetic-field strength, 

— B=(bb:Vu)B + r]AB-T 1 \Vb\ 2 B, (28) 
at 

and drop both the resistive terms and the second derivatives of 
the velocity field (bending terms), we might observe that 



Ck ,b - 



(26) 



-ln(B^ 1 / 2 ) = i(nn:Vu). 

dt v ' 2 V ; 



(29) 



The formal solution in the comoving frame is 

\n(BK l l 2 )(t)= l - J'dt'im:Vu(t')=^C n (t), 
while for the field-strength, we have 



lnfi(f): 



bb:Vu(t')=-( h (t). 



(30) 



(31) 



(K 2 ){B 2 ) ■ 

In all runs, it is found to be within 6% of its minimum possible 
value of -1 (Table 1). Furthermore, we notice that the mean 
square tension force [see eq. (1)] is 

(\B-VB\ 2 ) = (B 4 K 2 )+(B 4 \\7 ll B/B\ 2 ), (27) 

so the fact that K 2 m;i <~ (|B • VB| 2 /B 4 ) is large (and grows with 
Pr m ) while kj = (|B • VB| 2 )/(B 4 ) - k 2 v also indicates that the 
curvature in the regions of growing field remains comparable 
to the inverse scale of the flow. 8 

Let us give an illustration of the anticorrelation between B 
and K for a typical snapshot of the field. In Schekochihin et 
al. (2002b), we noted that this anticorrelation was manifest in 
cross sections of field-strength and curvature. Figure 5 does, 
indeed, show that field lines are straight (and direction revers- 
ing) in the areas of strong field. The same point can be made 
by scatter plots of B versus K during the kinematic stage of 
our run S4 (Fig. 8). Since the distributions of both magnetic 
field (§ 3.4) and curvature (§ 3.2.4) are intermittent, the scatter 

8 The second term in eq. (27) represents the contribution from the mirror force V11B/B, which is also large only in the bends, where its rms value grows at the same 
rate as K rms (Schekochihin et al. 2002b). The pdf of the mirror force has a power tail with the same scaling as the pdf of curvature (§ 3.2.4) and is, therefore, also 
dominated by the outer scales (A. A. Schekochihin 2003, unpublished). 

9 In Schekochihin et al. (2002c), we argued that BK ~ const based on a flux-conservation argument, which, however, involved some ad hoc assumptions about the 
fold geometry. Constantin et al. (1995) and Brandenburg et al. (1995) also argued in favor of BK ~ const on the basis of an alternative form of eq. (25). The numerical 
evidence presented here and in § 4.1 appears rather to support BK 1 / 2 = const. It is, however, possible that in the regions of strong field and low curvature, the relation 
between B and K is closer to B ~ l/K. 



It can be shown (A. A. Schekochihin 2002, unpublished), for 
the Kazantsev model velocity (7) and in neglect of bending and 
diffusion, that the joint pdf of 0,(0 A and („(t)/t is exactly the 
same as the joint pdf of C,\(t)/t and C2WA, the finite-time Lya- 
punov exponents corresponding to the stretching and the "null" 
direction of the flow (see § 3.3 for a quick overview of the rele- 
vant definitions). This means that, while Q,(t ) increases linearly 
with time and is responsible for the field stretching, ( n (t) fluc- 
tuates around zero. Equation (30) then suggests that BK 1 / 2 is 
a special combination in which the effect of stretching is can- 
celled. 9 Least-squares fits performed on log-scatter plots of B 
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versus K such as Fig. 8 give BK a = const with a ~ 0.47, which 
gives a measure of both the detailed anticorrelation between B 
and K and the extent to which it can be understood via essen- 
tially geometrical arguments such as the one we have just pre- 
sented. 10 Clearly, these should break down for curvatures close 
to either the inverse flow scale (K <~ k v ) or the inverse resistive 
scale (K <~ k v ), but it is nontrivial that they appear to work quite 
well away from the cutoffs. 
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FIG. 8. — Scatter plot of B vs. K at ; = 8 during the kinematic stage of our 
ran S4 (the 256 3 data were thinned out by a factor of 1000). 



3.2.4. The Curvature Distribution 

We have seen that the large values of K lms are due to the large 
values of curvature in the bends, where the field itself is weak. 
Figure 8 indicates that these only occupy a small fraction of the 
volume (in most places, fields are strong and straight). In order 
to ascertain that this is true, as well as to get a more detailed 
statistical description of the field-line geometry, Schekochihin 
et al. (2002b) found the pdf of curvature analytically: 



Pk(K) = - 



K 



7 [l + C^/X) 2 ] 10 / 7 ' 



(32) 



where K* = (2ka/1 m) 1 ! 2 [see eq. (8) for definitions of K2 
and k 4 ]. This solution means that in most of the volume, 
the field-line curvature is comparable to the inverse eddy size 
(K <~ <~ k v ), while the distribution of curvatures in the bends 
is characterized by a power law P K (K) ~ K~ 13 ^. This scaling 
is reproduced very well in our simulations (see Fig. 9). Note 
that we find all curvature-related quantities to be quite well 
converged already after a relatively short running time (unlike 
the quantities containing field-strength, which fluctuate very 
strongly). We believe that this is due to statistics of curvature 



being a geometrical property of the field lines independent of 
the fluctuating stretching rates (see footnote 10). 

A stationary power-like pdf of K is possible because curva- 
ture, unlike field strength, has an explicit dependence not only 
on Vu but also on VVu. The second derivatives of u enter as 
a source term in eq. (25) and are responsible for bending the 
fields. Thus, the flow scale is explicitly present in the curvature 
equation and, consequently, in the curvature statistics. 

Formally speaking, integer moments of the distribution (32) 
diverge: in the diffusion-free calculation of Schekochihin et al. 
(2002b), while the pdf converges to the stationary profile (32), 
all moments (K n ) grow exponentially without bound. In the 
problem with diffusion, the power tail of the curvature pdf is 
cut off at the resistive scale (K <~ k v ) and the moments saturate 
at values that scale with Pr m : e.g., ^ rms ~ fcyPr^ 7 asymptoti- 
cally with Pr m — > oo. The numerical results demonstrate the 
(nontrivial) fact that the curvature statistics above the resistive 
scale are not affected by the presence of diffusion. In other 
words, there is no resistive anomaly for the curvature: r\ — > +0 
and r\ = give the same result. Note that the use of hyperdif- 
fusion does not change the scaling of Pk(K) (Schekochihin et 
al. 2002c) — another indication that curvature statistics do not 
depend on the dissipation mechanism. It can be seen in § 3.4 
that, unlike the curvature, the field strength has statistics that 
are crucially influenced by the resistive cutoff. 



PDFs of curvature 

kin. diffusive 

saturated 



K-13/7 




K/2TT 

FIG. 9. — Curvature pdf's for run S4 during the kinematic regime (with diffu- 
sion) and in the saturated state. Curvature pdf's for all other runs are similar. 

A note of caution is in order. We have found in the course 
of our numerical investigations that the geometric properties of 
the dynamo-generated field, while very well behaved in terms 
of convergence in time, are quite sensitive to spatial resolu- 
tion. Specifically, with our spectral code, if the magnetic dif- 



10 The underlying geometrical nature of the curvature-related statistics becomes especially clear in 2D. The curvature and field-strength in 2D are related by 
BK l l 3 = const (again neglecting diffusion and bending). This follows from the fact that, in 2D, n _L b, so, b being the stretching direction, n must necessarily 
be the compressive one: incompressibility then requires f„ = whence dln(B K 1 / 3 )/dt = 0. A statistical calculation for the Kazantsev model shows that the pdf of 
\n(BK [ / 3 ) is a <5 function (A. A. Schekochihin 2002, unpublished). A purely geometrical consideration of the field lines also gives BK 1 / 3 = const (Thiffeault 2004). 
Furthermore, it turns out that the power tail of the curvature pdf in 2D derived by Schekochihin et al. (2002b) for the Kazantsev dynamo model, Pk(K) ~ K~ 5 I 3 , can 
be reproduced for the curvature distribution along a generic parabola (J.-L. Thiffeault 2002, private communication). It is an open question whether similar purely 
geometric treatment is possible in 3D. 
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fusivity 1] is not large enough to ensure a fairly large separa- 
tion between the resistive scale and the dealiasing cutoff, the 
folded structure is polluted by spurious small-scale "ringing" 
in the field. The method of computing curvature is also impor- 
tant. Better-quality results are obtained by using the formula 
K = (B VB) ■ (l-bb) /B 2 , than by directly differentiating the 
unit vector (K = b Vb). 

3.3. How Is the Small-Scale Dynamo Possible? 

The folded-structure diagnostics considered in § 3.2 tell us 
what field configurations are generated by random stretching. 
Diffusion did not figure prominently in the theoretical discus- 
sion that we have offered in support of our numerical results. 
Its only effect has been to set the minimum scale of field varia- 
tion (reversals and curvature). However, if the fields reverse at 
the resistive scale, why are they not destroyed by diffusion? In 
other words, how is the small-scale dynamo possible at all? 

Let us consider the simplest imaginable model of turbulent 
stretching: a velocity field random in time and exactly linear in 
space (Zel'dovich et al. 1984). In this model, any physics that 
depends on the finiteness of the flow scale is lost: for example, 
the stationary curvature pdf (§3.2.4) cannot be derived for an 
exactly linear velocity (a more serious shortcoming will be the 
incorrect result for the growth rates of (B n ): see §3.4). The 
fields produced by a linear flow still have a folded structure, but 
there is no bending, so the folds can be arbitrarily long. 

If the velocity field is linear, eq. (2) is solved (in the comov- 
ing frame) by the Ansatz 



B(f,x) = 



(Hcq 



■B(f,k )e'' k( '' ko)x , 



(33) 



where k(0,ko) = ko, and B(0,ko) = B (ko) is the Fourier trans- 
form of the initial field. After direct substitution into eq. (2), 
we find the following solution 



dx' 

B '( f > k o)=7r^( f )Bo(ko)exp 

ox 



7] / <if'A: 2 (f,ko) 



dxT 



(34) 



(35) 



where Xrj — 
the flow, i.e 



x(f ) is the transformation of variables induced by 
x(f) = xo(f) + Jq dt 'u(f ',x(f '))■ For a linear velocity 
field, dx' /dx" } ' depends on time only. 

Solution (34) is a generalization of the Cauchy solution of 
the ideal induction equation. It can be used to express the statis- 
tics of the magnetic field in terms of its initial statistics and the 
statistics of the (contravariant) metric associated with the flow 



g' J (t) = 



dx 1 dx-* 
d^dx% 



(36) 



The metric contains all the statistical information necessary to 
describe random advection of arbitrary tensor fields by the am- 
bient velocity field (e.g. Boldyrev & Schekochihin 2000). Since 
g'j is symmetric and positive definite, it has a set of real positive 
eigenvalues {e^ 1 , <?^ 2 , e'*- 1 } and an orthonormal basis {ei,e2,e3}. 
This is the Lyapunov basis of the flow, and Q{t)/t are the finite- 
time Lyapunov exponents. The quantities can be considered 
ordered so that d > £2 > (3- Incompressibility requires that 
Ci + C2 + C3 = at a U times. Therefore, Ci > corresponds to the 
stretching direction, (3 < to the compressive direction, and 

" In other words, fi = fj, where £(, is defined in eq. (31). As we noted in §3.2.3. 
direction and £2 = C« [ e 1- (30)] for a linear velocity field. 



£2 to the "null" direction. The exponent (2 can be either posi- 
tive or negative. In time-reversible flows (e.g., in the Kazantsev 
model), (£2) = (in real turbulent flows, it is usually positive). 
Note that Lyapunov eigenvectors stabilize exponentially fast in 
time, while the finite-time Lyapunov exponents Q/t converge 
to constant values as ~ l/t (Goldhirsch et al. 1987). 

Equation (34) implies that B aligns with the stretching direc- 
tion ei . In an ideal fluid, 

3 2 



fi z ~exp(G) 



(37) 



for a typical realization. 11 While B transforms as a vector 
[times the exponential factor due to diffusion: see eq. (34)], 
k transforms as a co vector [eq. (35)]. Therefore, it aligns with 
the stretching direction of the covariant metric, which is the 
inverse of g'K Its stretching direction is then the same as the 
compressive direction of g l ->, with exponent -(3 > 0, so k aligns 
with the compressive direction £3 and grows exponentially in 
time. Therefore, most modes B(f,ko) [eq. (34)] decay super- 
exponentially. The only modes that are not thus suppressed 
are those for which the angle between the initial wavenum- 
ber ko and the compressive direction is very close to 90°, with 
the window of allowed angles narrowing exponentially fast in 
time (Zel'dovich et al. 1984). Since klB and B || ei, we 
have k || 62, i.e., the surviving fields can have reversals only 
in the "null" direction. In contrast, in two dimensions, the field 
aligned with the stretching direction must necessarily reverse in 
the compressive direction, so the folds are destroyed by diffu- 
sion (see simulations by Kinney et al. 2000). These statements 
(illustrated by Fig. 10) explain qualitatively why small-scale 
dynamo works in three dimensions but not in two dimensions 
(Zel'dovich 1957). Even in three dimensions, diffusion leads 
to suppression of most modes. Zel'dovich et al. (1984) showed 
that the remaining ones, with the correct alignment of k and B, 
are sufficient to make the total magnetic energy grow as 

exp[(6 00 2|. 



B 1 



(38) 




2D 



Compression. 
S *3 



*" e-, 



Stretching 



Comp -ession 




FIG. 10. — Fold alignment in three and two dimensions. 

In the linear-velocity model, a system of infinite size can 
contain an infinite number of initial wavenumbers ko (no new 
modes are produced by a linear velocity field). Only an expo- 
nentially small fraction of ko values contribute to the growing 
field at any given time. In a finite system, only a finite num- 
ber of modes can exist, limited by the system-size and resistive 
cutoffs. In a linear velocity field, they would all eventually be 

it is also possible to show that the curvature K = b ■ Vb aligns with the "null" 
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suppressed. In a finite-scale flow, besides stretching, there is 
also bending of folds, which effectively means that new modes 
are made. Thus, the effect of the finiteness of the flow scale is 
that modes are continually resupplied. 

One might wonder how the Zel'dovich et al. (1984) small- 
scale dynamo mechanism relates to the classic "stretch- 
twist-fold" (STF) picture of the fast dynamo (Vainshtein & 
Zel'dovich 1972; Childress & Gilbert 1995). The mechanism 
outlined above makes the small-scale dynamo possible as a sta- 
tistical effect: the random flow produces certain typical struc- 
tures (folds), then diffusion selects a subset of these structures 
that turns out to be sufficiently large on the average to give rise 
to energy growth. No specific form of the flow is prescribed: 
three-dimensionality, randomness in time, and spatial smooth- 
ness are all that is required. In contrast, the STF picture shows 
what might happen to a flux loop in three dimensions if a par- 
ticular set of transformations is applied to it. Vainshtein et al. 
(1996) studied numerically the evolution of a single circular 
field line in a chaotic flow that was designed to favor STF trans- 
formations. They observed that only at the very beginning of 
the evolution could the line be seen to actually stretch, twist, 
and fold. It then quickly became randomly tangled and, while 
the field did grow exponentially, it was not possible to claim 
that the mechanism of this growth was successive stretching, 
twisting, and folding. Note also that in order to twist the field 
lines, the flow must possess pointwise helicity (net helicity is 
not required because the sense of the twist can be arbitrary). 
However, numerical experiments by Hughes et al. (1996) sug- 
gest that flows with zero pointwise helicity can still generate 
small-scale fields. Thus, while the STF mechanism is, in prin- 
ciple, a viable dynamo, it remains unclear whether a generic tur- 
bulent flow contains enough STF realizations to amplify mag- 
netic energy this way, whether this is the statistically prevalent 
kind of dynamo, and what its relation to the Zel'dovich et al. 
(1984) mechanism is. 

It would be worthwhile to make careful measurements of the 
fold alignment with respect to the local Lyapunov basis in nu- 
merical simulations. We leave such a study outside the scope 
of this paper. 

3.4. Self-Similar Growth and Intermittency 

The subviscous range is bounded on the IR side by the vis- 
cous scale k~ l (the flow scale) and on the UV side by the re- 
sistive scale k' 1 . Taking k v — ► amounts to assuming a linear 
velocity field; k n — > oo is the diffusion-free limit. Both limits 
are singular for at least some of the statistical quantities of in- 
terest. For the magnetic-energy spectrum, the resistive cutoff k n 
is explicitly present, but the IR cutoff k* ~ k v only enters via 
a small correction to the growth rate [eq. (12)]. On the other 
hand, the flow scale is present explicitly in the quantities re- 
lated to the field structure: thus, e.g., the field-line curvature 
is allowed to have a stationary pdf only because of the broken 
scale invariance on the IR side [K* ~ k v , eq. (32)], while diffu- 
sion simply cuts off the pdf at large K. In both cases, analytical 
theory is facilitated because one of the cutoffs does not matter. 

For the field-strength distribution, both cutoffs turn out to be 
important. Analytical theory exists only for the diffusion-free 
case (k v = oo, § 3.4. 1) and for a purely linear velocity field with 
diffusion (k v = 0, k v finite, § 3.4.2). Neither of the two limits 
correctly describes intermittency of the magnetic field gener- 
ated by a finite-scale flow (§ 3.4.3). The linear-velocity model 
(reviewed in § 3.3) is very helpful in understanding the way the 
small-scale dynamo works in the presence of diffusion. The in- 



adequacy of this model as a theory of intermittency has not until 
now been properly appreciated, so we would like to emphasize 
the importance of the results reported in § 3.4.3. 

3.4. 1 . The Diffusion-Free Case 

The situation in the diffusion-free regime is straightforward. 
Equation (2) with r\ = admits the formal solution (31) forB(f) 
in the comoving frame, so, on the basis of the Central Limit 
Theorem, B should have a lognormal distribution. As usual, 
the Kazantsev model (7) allows for an exact solution. The re- 
sulting pdf of B is, indeed, lognormal, with a growing mean 
(ln[B(f)]) = (3/4) K 2 t and dispersion D = K 2 t [see eq. (41)]. As 
the width of this lognormal profile grows with time, the field 
becomes increasingly more intermittent: the moments of the 
field strength are 

(B n ) ocexp[(l/4)K 2 n(n + 3)f]. (39) 
This formula can be obtained both from direct averaging of the 
induction equation (e.g., Schekochihin & Kulsrud 2001) and 
from averaging eq. (37) over the statistics of the finite-time Lya- 
punov exponent (i(t)/t (Chertkov et al. 1999). 

A good conventional measure of intermittency is the kurtosis 

{B 4 )/{B 2 } oc exp(2 K2 f), which can be roughly interpreted as 
the inverse volume-filling fraction. The kurtosis and all other 
normalized moments (B mn )/(B m ) n grow exponentially. This 
means that the dominant contribution to each moment (B n ) 
comes from a different substructure, which occupies an ex- 
ponentially decreasing fraction of the volume compared with 
those contributing to lower moments. This is because the fields 
everywhere are stretched exponentially, but with a fluctuating 
stretching rate C,\(t)/t [eq. (37)], so any occasional difference 
between substructures tends to be amplified exponentially. 

3.4.2. A Model with Diffusion and a Linear Velocity Field 

The situation changes dramatically when diffusion is intro- 
duced. Exact theory for one-point statistics in this case is prob- 
lematic even for the Kazantsev velocity, because of the clo- 
sure problem associated with the diffusion term. The prob- 
lem is solvable, however, for the particular case when the ve- 
locity field is exactly linear and the system domain is infinite 
(Zel'dovich et al. 1984; Chertkov et al. 1999; Nazarenko et al. 
2003; West et al. 2004). The general outline of the solution was 
given in § 3.3 leading to the growth of a typical field realization 
according to eq. (38). Chertkov et al. (1999) worked out the 
moments (B 2n ) averaged over the (Gaussian) distribution of Q 
for the Kazantsev model and found 

(B 2n ) cxPr^ /4 exp[(3/16)K 2 n(« + 4)r]. (40) 
Thus, it appeared that, in the diffusive regime, while the spe- 
cific expression for the moments changed, the intermittency 
continued to increase exponentially in time just as it did in the 
diffusion-free case. 

3.4.3. The Finite-Scale Flow: Self-Similarity 

We already saw in § 3.2 that the linear- velocity model is 
not sufficient to describe the field structure. The intermittency 
properties of the field strength are also affected — and, indeed, 
drastically changed, — by the presence of the finite system 
scale. 

In an infinite system, intermittency can grow with time be- 
cause for ever higher moments (B n ), ever smaller sets of sub- 
structures can always be found in which the field has exponen- 
tially outgrown the rest of the system and that, therefore, domi- 
nantly contribute to (B n ). By contrast, in a finite system, only a 



SIMULATIONS OF THE SMALL-SCALE TURBULENT DYNAMO 



13 



(a) 



i 1 i 1 i 1 r~ 

Run S4: Pr =500 



'! 



I! 



(B 4 > 1 / 2 /(B 2 ) 
(B 6 ) 1 / 3 /(B 2 ) - 
<u 4 ) 1 / 2 /(u a ) _ 



1 ' I' I, 

1 \/V 




kinematic (diffusive) nonlinear 



z4 79.d 

70 



(b) 

° kinematic (ciffusive) 
■ saturated 



i 

Runs S0-S6 



Pr o. 



i 1 i 



/Kiii d 
z476-d 
z479.d 
z584.d z471.d 
z585.d z472.d 
z586.d z634.d 
z587.d 



10 3 



Pr 



FIG. 11. — (a) Evolution of normalized moments of the magnetic-field strength for run S4. The (square root of) kurtosis of the velocity field is given for compari- 
son, (b) Kurtosis (B 4 ) /{B 2 ) 2 of the magnetic-field strength during the kinematic stage of runs S0-S4 and in the saturated state of runs S1-S6. The values plotted are 
listed in Table 1. 



finite number of exponentially growing substructures can exist, 
so the contribution to all moments must eventually come from 
the same fastest-growing one. The statistics of the field should, 
therefore, be self-similar, with (B n ) growing at rates propor- 
tional to n, not n 2 , and all normalized moments (B mn ) / (B m ) (B n ) 
saturating. 12 

In our simulations, there is unambiguous evidence for the 
self-similar evolution. After initial diffusion-free growth, the 
normalized moments saturate (Fig. 11a). The pdf of the 
magnetic-field strength becomes self-similar: namely, the pdf 's 
of B/B rms collapse onto a single stationary profile throughout 
the kinematic stage of the dynamo (Fig. 12). The large-B tail 
of the pdf of B/B lms is reasonably well fitted by a lognormal 
distribution. Specifically, suppose that the pdf of z = InB is 

-1/2, r / „ „>\2 



P z (z) = (nDT 



(41) 



2 exp[-(z-(lnB»>]. 

Then (B n ) ocexp [(In B)n + Dn 2 / 4], so D = ln((B 4 ) l/2 / (B 2 )). In 
the diffusive regime, D = const and the field-strength statistics 
become self-similar: the pdf of ( = ln(B/B rms ) is stationary, 

P C (C) = (7rD)- 1 / 2 exp[-(C + Z)/2) 2 /£'] • (42) 
The lognormal fit in Figure 13 is obtained by calculating D 
from the numerical data and comparing the profile (42) with 
the numerically calculated pdf. The fit is qualitative but decent, 
considering the simplicity of the chosen profile (42) and large 

statistical errors in determining (B 4 ) / (B 2 ) (Fig. lib). 

The pdf at small B appears to be power-like, BP(B) ~ B 2 - 7 - 2 - 8 
in all regimes (Fig. 12b). This tail extends to the limit of 
our diagnostic's resolution [BP(B) ~ 10~ 5 at B ~ 10" 3 B rms , not 
shown], but there may be an unresolved lognormal tail at even 
smaller B. It is clear that, while the large-B tail describes the 
straight segments of the folds, the small-B tail gives the field- 
strength distribution for the weak fields in the bends. 

12 This effect was not properly appreciated in our discussion of the kinematic dynamo in Schekochihin et al. (2002c). Note that for the problem of passive-scalar 
decay, self-similar behavior was found in 2D maps (the so-called "strange mode": see Pierrehumbert 1994, 2000; Sukhatme & Pierrehumbert 2002) and even in 
scalar-mixing experiments (Rothstein et al. 1999). Numerical simulations of the advection-diffusion equation in two dimensions have shown the same property 
(Fereday & Haynes 2004). In map dynamos studied by Ott and coworkers (Ott 1998), all moments of the field also grew at the same rate. 



Note that, in space, intermittency of the field strength means 
that the growing fields do not uniformly fill the volume. More 
specifically, intermittency is often associated with the presence 
of "coherent structures" that have disparate spatial dimensions 
For the dynamo-generated magnetic fields, they are the folds, 
for which the disparate dimensions are their length and the 
field-reversal scale. A linear velocity field allows the folds to be 
elongated indefinitely, thus giving rise to arbitrarily large aspect 
ratios between the two scales. In reality, the length of the folds 
cannot be larger than the scale of the flow because of the bend- 
ing of the folds. The aspect ratio is then bounded from above 
by the maximum allowed scale separation ~ Pr^ 2 . 

The hypothetical lognormal pdf (41) is self-similar only if 
its dispersion D does not depend on time. In contrast, in the 
diffusion-free regime, we had D ~ 7?, where 7 ~ Vu is the 
stretching rate. In the case of rj > and linear velocity field, the 
formula for (B n ) derived by Chertkov et al. (1999) [eq. (40)] is 
also consistent with a lognormal distribution for which D <~ 7?. 
Both results are only valid transiently, during the time that it 
takes magnetic fluctuations to reach the resistive scale (in the 
former case) or the system (flow) scale (in the latter case). 
Since the scale separation is ~ Pr^ 2 and the spreading over 
scales proceeds exponentially fast at the rate <~ 7 (§ 3.1), the 
time during which intermittency increases is f* ~ 7 _1 lnPrJ/ 2 . 
Physically, this is the time necessary to form a typical fold 
with length of the order of the flow scale and field reversals 
at the resistive scale — starting from either a flow-scale or a 
resistive-scale fluctuation. We might conjecture that this time 
determines the characteristic magnitude of the dispersion in the 
self-similar regime: D ~ 7?* ~ lnPr^ 2 , which implies that the 

kurtosis (B 4 )/(B 2 ) = exp(2D) increases with Pr m in a power- 
like fashion. The specific power law depends on prefactors that 
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FIG. 12. — (a) Evolution of the pdf of the magnetic-field strength for run S4. The pdf's for the diffusion-free regime (0 < t < 2) are given at time intervals 
of Af = 0.4. The subsequent evolution (2 < t < 40) is represented by pdf's at time intervals At = 4. (b) Collapse of the field-strength pdf's in the kinematic stage 
of ran S4 onto a self-similar profile. These are pdf's of B/B rms . The pdf in the saturated state is shown for comparison. 



may be nonuniversal. The same holds for other normalized mo- 
ments. Figure 1 lb is an attempt to test this hypothesis for a 
sequence of simulations with increasing Pr m . The fluctuations 
in the kinematic regime are very large (see Fig. 11a), so the er- 
ror bars are too wide to allow us to claim definite confirmation 
of the power-like behavior, but our results are consistent with a 

A caveat is in order. At our resolutions, the lognormal fit 
is, in fact, not the only one that is compatible with the numer- 
ical evidence. Stretched-exponential and even steep-power-tail 
[P(B) ~ ZT 4 - 5 ] fits of comparable quality can be achieved. Note, 
however, that a stretched exponential would not be compatible 
with a Pr m -dependentkurtosis. 

We emphasize that the self-similarity reported here is statisti- 
cal, not exact. Namely, it does not imply that the magnetic field 
is simply a growing eigenmode of the induction equation (2). 
Such an eigenmode does exist for some finite-scale nonrandom 
(and time-independent) flows and maps (Childress & Gilbert 
1995). The self-similar pdf we have found here is a natural 
generalization of this eigenmode dynamo to random flows. 

A testimony to the statistical nature of the self-similarity is 
the tendency of the growing field configuration to undergo sud- 
den disruptions when the self-similar profiles are rapidly de- 
stroyed and then regenerated by the dynamo. These disruptions 
are manifested, e.g., in sudden destruction of the large-fe tails 
of magnetic spectra (Fig. 4a). These events are correlated with 
large downward fluctuations of the velocity field. Indeed, when 
such a fluctuation occurs, the stretching rate 7 <~ Vu drops, 
the resistive scale l n ~ (?y/7) 1//2 increases, and the whole field 
structure with direction reversals at the previously established 
smaller resistive scale is quickly dissipated. Since it is the 
stronger straight direction-reversing fields that are destroyed, 
the tail of the field-strength pdf is suppressed. The self-similar 
statistics are quickly reestablished as the magnetic field adjusts 
to the new stretching rate. Of course, the disruptions are noth- 
ing more than extreme events in what is a continuous series of 



fluctuations. Their presence shows the degree of disorder be- 
hind the self-similar statistics. 




B/B n 



FIG. 13.— ThepdfofB/B rms for run S4 (averaged over the kinematic dif- 
fusive stage) and the lognormal profile (42) with the same D. Also given is the 
pdf in the saturated state. The pdf's for all other runs are similar. 

3.5. Summary of the Kinematic Dynamo 

We have established the following properties of the kine- 
matic regime. 

1 . The bulk of the exponentially growing magnetic energy 
is concentrated at the resistive scale l v <~ Pr~ 1,/2 Z, y . The 
magnetic-energy spectrum peaks at k n <~ and in the 
subviscous range k v <C k <C k v , it has a positive spectral 
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exponent consistent with the theoretical prediction k 3 / 2 
(Kazantsev 1968; Kulsrud & Anderson 1992). 

2. The fields generated by the dynamo have folded struc- 
ture: they reverse direction at the resistive scale with the 
field lines remaining straight (and fairly well aligned) up 
to the scale of the flow. It is the direction reversals that are 
responsible for the magnetic-energy spectrum peaking at 
the resistive scale. In the regions of curved field (bends), 
the field is weak. The field's strength and its curvature 
are, thus, anticorrelated. The theory of Schekochihin et 
al. (2002b) is broadly confirmed. 

3. The growing direction-reversing fields are intermittent in 
space. The field-strength distribution is consistent with 
a lognormal profile and evolves self-similarly with time. 
The volume fraction occupied by the growing fields de- 
creases with increasing Pr m in a power-like fashion. 

The Zel'dovich et al. (1984) mechanism of fold alignment 
(§ 3.3) explains the existence of the dynamo. 

4. THE NONLINEAR SATURATED STATE 

The kinematic growth of the magnetic energy must culminate 
in a saturated nonlinear state. When does the nonlinearity [the 
Lorentz tension force in eq. (1)] become important? We showed 
in §3.2 that the kinematic dynamo produced folded fields for 
which B • VB ~ k v B 2 . This means that the tension force can 
start to oppose the purely hydrodynamic terms in eq. (1) once 
B • VB ~ u • Vu B 2 ~ u 2 v , i.e., when the magnetic energy 
becomes comparable to the energy of the viscous-scale eddies. 
Note that this is a nontrivial statement that only holds for the 
folded field: in an unstructured field with the same k rms ~^,we 
would have had B • VB ~ k v B 2 and the nonlinearity would have 
set in at much lower field energies B 2 ju 2 v ~ k v /k v ~ Pr" 1 ^ 2 . 
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Fig. 14. — Saturated values of the kinetic and magnetic energies vs. Pr m 
values plotted are listed in Table 1 . 
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equal kinetic and magnetic energies. We find that the satu- 
rated values of (B 2 ) grow with Pr m toward an apparent constant 
asymptotic value somewhat above (u 2 ) (Fig. 14, cf. Brummell 
etal. 2001). 

Although a number of heuristic models of back-reaction have 
recently been proposed (Subramanian 1999, 2003; Kim 1999, 
2000; Boldyrev 2001; Nazarenko et al. 2001; Schekochihin et 
al. 2002c,d, 2004b, and § 4.3), the exact mechanism of the non- 
linear saturation in isotropic MHD turbulence is not currently 
well understood even on the qualitative level. One must, there- 
fore, probe the nonlinear physics by careful numerical analysis. 
As always in turbulence problems, only statistical statements 
can be expected to have any degree of universal applicability. 
The statistical averages are replaced by time averages, so, in 
order to collect statistics, simulations must be run for a long 
time after the statistical steady state has been achieved. We 
present statistical results averaged over 20 time units, after all 
discernible systematic change of the nonlinear state ceased. 

We first establish the survival of the folded field structure in 
the saturated state (§ 4.1), then look at the magnetically induced 
subviscous velocity fluctuations (§ 4.2), and, finally, study the 
magnetic-energy spectrum, for which a nonlinear saturation 
model is described and tested (§ 4.3). Intermittency in the satu- 
rated state is considered in § 4.4. A summary is given in § 4.5. 



FIG. 15. — Field structure at t - 
nonlinear counterpart of Fig. 5. 




76 in the saturated state of run S4. This is the 



In runs S1-S6, the viscous scale is the same as the forcing 
scale, so the above arguments imply saturation with roughly 



4. 1 . Persistence of the Folded Structure 

A key fact that we learn from numerical simulations is that 
the field stays folded after its energy saturates (Fig. 15). All 
main features of the folded structure observed in the kinematic 
regime (§ 3.2) carry over to the saturated state: 

1 . Field varies along itself at the scale of the flow. Its char- 
acteristic parallel wavenumber [eq. (13)] is k\\ ~ k v and 
does not depend on Pr m , while &bxj and £ lms [eqs. (14- 
15)] increase with Pr m (Fig. 16, cf. Brandenburg et al. 
1995; Brummell et al. 2001). The scaling k rmti ~ k^xj ~ 
Pr,!/ 2 appears to be recovered asymptotically with Pr m , 
but is much less well satisfied at our finite Pr m than it 
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was in the kinematic case (cf. Fig. 6b). This is consis- 
tent with the flattening of the magnetic-energy spectrum 
discussed in § 4.3. 

2. The magnetic-field strength and curvature remain anti- 
correlated with the correlation coefficient C k ,b [eq. (26)] 
within 4% of -1 (Table 1) and the log-scatter plots of 
B versus K least squares fitted by BK a = const with 
a ~ 0.45 (Fig. 17; cf. § 3.2.3). 

3. The curvature pdf in the saturated state has a power tail 
with scaling very close to the kinematic K~ 13 / 7 (the scal- 
ing in the saturated case may be slightly steeper than in 
the kinematic case: see Fig. 9). The bulk of the pdf is 
concentrated at the flow scales (cf. § 3.2.4). 
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FIG. 16. — Averaged values of various characteristic wavenumbers (denned 
in §3.2) vs. Pr m in the saturated state of runs S1-S6. The values plotted are 
listed in Table 1. Note that the value of for run S6 is too large as a result of 
this run being somewhat underresolved. 

Physically, it is clear that the preservation of the folded struc- 
ture is inevitable. The large viscosity of the medium does not 
allow the setting up of the kind of detailed small-scale flows 
that would be required to unfold the field without destroying it. 
Furthermore, a wholesale evacuation of the small-scale degrees 
of freedom would lead to an unstable situation in which any 
significant velocity fluctuation would set in motion small-scale 
field generation and refolding. The latter point was checked 
numerically by Maron et al. (2004). 

The persistence of the kinematic scaling of the tail of the cur- 
vature pdf also makes sense. The tail describes the curvature 
distribution in the bends of the folds, where the magnetic field 
remains weak even when field energies in the straight segments 
of the folds approach equipartition levels. Therefore, the situ- 
ation in the bends is quasikinematic, so the curvature statistics 
are almost the same as in the kinematic limit. 

We note one feature of the field structure in the saturated 
state that differs from the kinematic regime: k#.j [eq. (17)] in- 
creases somewhat (Fig. 6a) and appears to scale with Pr m in the 
same way as k ms and kn x y. we find k rms : &bxj : &b j ^5:3:1 
(Fig. 16). This indicates that the direction-reversing fields are 



less well aligned than in the kinematic case. Geometrically, 
one can interpret this effect as a transition from flux sheets to 
flux ribbons. While intuitively such a transition might be at- 
tributable to increased mixing efficiency of the velocity field 
in the saturated state (see §4.3.1), we do not currently have a 
satisfactory theoretical picture of this nonlinear effect. 

Saturation requires some form of suppression of stretching 
by the tension force in eq. (1). Since the stretching motions 
are large-scale, the back-reaction must have a large-scale co- 
herence. The folded fields, while formally small-scale, possess 
such a large-scale coherence in the form of slow (in space) par- 
allel variation. The persistence of the folded structure in the 
saturated state is crucial for the saturation model in § 4.3. 
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FIG. 17. — Scatter plot of B vs. K at t = 76 in the saturated state of our run S4 
(the 25 6 3 data were thinned out by a factor of 1000). 



4.2. The Tension Force and the Subviscous Velocity Spectrum 

The first traces of nonlinearity appear already in the kine- 
matic stage in the form of the small-scale modification of the 
velocity spectrum. For k 3> k v , the Navier-Stokes equation (1) 
reduces to the following force balance (in k space): 



Mt 2 u(k)~-/k/?(k) + F(k), 



(43) 



where F = B • VB is the tension force. Pressure is determined 
from the incompressibility constraint, k • u(k) = 0, and we get 
for the velocity spectrum at k 3> k v 



E(k) = iy dQk^dutk)! 2 ) 



-^/^( F ( k )-( S -F)- F( - k) )- (44) 

By the same argument as in §3.2.2, the angle integral on the 
right-hand side is ~ k°, so E(k) ~ k~ 4 '. Figure 18 shows that 
the velocity spectrum does indeed develop such a tail at sub- 
viscous scales already during the kinematic-dynamo stage. The 
tail thickens as the tension force grows (F oc B 2 ). 
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FIG. 18. — Evolution of the kinetic-energy spectrum for run S4. We plot the 
spectrum in the purely hydrodynamic regime (B = 0), the spectrum in the satu- 
rated state, and the evolving spectra for 16 < t < 40 at time intervals of At = 2. 
The corresponding evolution of the tension spectra is given in Figure 7a. 

In the saturated state, the tension spectrum remains flat, as 
does M 4 (k), and the relation T(k) ~ k^M^k) persists: indeed, 
since the folded structure is preserved, the argument given in 
§ 3.2.2 should continue to apply. The velocity spectrum keeps 
its k~ A tail, and we find that 



E{k)> 



2 1 



3 v 2 k 4 

is quite well satisfied (Fig. 19). 
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FIG. 19. — Tension spectrum for run S5 (saturated state). For comparison, 
we also plot the kinetic-energy spectrum E(k) compensated by (3/2)u 2 k 4 and 
the spectrum of B 2 , multiplied by k 2 [as defined by eq. (13)]. The results for 
runs S1-S4 and S6 are analogous. 



Thus, the subviscous k~ A tail of the velocity spectrum is a 
passive feature. It contains very little energy and has nothing 
to do with the large-scale back-reaction mechanism that causes 
saturation. It is simply another signature of the folded structure 
and is present for both kinematic and nonlinear dynamo. 

4.3. Local Anisotropy and Steady-State Spectra 

Let us now turn to the key issue of the dynamo physics: 
how does the dynamo saturate? There must exist a nonlinear 
feedback mechanism that (1) maintains a statistically steady 
magnetic-energy level and (2) is compatible with the strongly 
fluctuating velocity field continuously stirred by the random 
forcing. Here we outline a model for the saturated spectrum 
we proposed recently (see Schekochihin et al. 2004b, where it 
is presented in a somewhat different form) and compare its pre- 
dictions with our simulations. 



4.3.1. Dynamo Saturation via Local Anisotropization 

Saturation can be achieved if stretching motions are nonlin- 
early suppressed. This does not have to mean complete elimina- 
tion of all turbulence: only the bb : Vu component of the rate- 
of-strain tensor leads to work being done against the Lorentz 
force and must, therefore, be suppressed. It is then natural that 
the velocity field should become anisotropic with respect to the 
local direction of the magnetic field. Since the magnetic field 
has a folded structure (§ 4. 1), such a local direction is defined by 
the tensor bb J ', which varies at the scale of the flow (direction 
reversals cancel). The back-reaction term B • VB in eq. (1) is 
quadratic in B, so it makes sense that the nonlinearly modified 
velocity should be indifferent to field reversals but anisotropic 
with respect to the direction of the folds. 

In order to model such an anisotropization, we build on the 
existing kinematic theory based on the Kazantsev velocity (7). 
Thus far, isotropic correlators K ,J (y) have been assumed. Now 
we let k'-> depend on the local direction of the folds. In the pres- 
ence of one preferred direction defined by the tensor b'b\ the 
correlator of an incompressible velocity field in k space has the 
following general form 

K ij (k) = K & (k, \/j,\)(6 ij -kkj) +K ( *\k, \n\)(b r b j 

+ n 2 kikj - fib'kj - Likib j ) , (46) 

where k = k/k, /j, = k • b. It is then possible (Schekochihin et al. 
2004b) to derive an equation for the magnetic-energy spectrum 
at k > k v that takes into account the dependence of the velocity 
statistics on the local direction of the field: 



dtM= ^dk 



(1 + 2a«)k 2 ^- -(1 +4<t_l + l0au)kM 

II Qfc II 



+ 2(cr_L + cr||)7 ± M-2r/fc 2 M, 



where 



7-L 



cr_L = 



d s k 



■fclKi_(k), 



1 
1 



Ull = 



d 3 k 
(2tt)3 
d 3 k 



7_L J (2iry 



kj K± (k) 7 
fcjKii(k), 



(47) 

(48) 
(49) 
(50) 
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and we have denoted k± = k(l -/i 2 ) 1 / 2 , ku = kfj,, 
K±(k) = l(S ij -b i b i )K ij (k) 



= \ [(l + ^^lMD + ^d-^ 2 )^^,^!)] ,(51) 



K||(k): 



2 

ifcW 7 (k) 



1 



(l-fj, 1 ) \K W (kM) + (l-»)K W (k,\» 



(52) 



In the isotropic case, k w = k (1) (£), K (a) = 0, which gives = 
2/3, (J|| = 1 /6, and -f± = (6/5)7 so e< J- (47) reduces to the stan- 
dard kinematic case [eq. (9)]. 

The derivation of eq. (47) depends on the assumption that the 
tensor b'b-' changes at the scale of the flow. This is only true in 
the regions where the fields are strong and straight. The regions 
of curved fields (bends of the folds) are not properly described 
by this formalism, but the fields there are weak (§ 4.1) and oc- 
cupy a small fraction of the volume compared to the strong 
straight ones (because of the steep power tail of the curvature 
pdf, Fig. 9). Since it is the stronger fields that dominate the 
back-reaction, we believe that our description is a reasonable 
one. 

With a zero-flux boundary condition imposed at some k* <~ 
k v and in the limit rj — > +0, eq. (47) has the eigenfunction 



M(k) = (const) e^k s Ko(k/kr,), 
(l+2a||) 7i l 1/2 



7 = 7± 



s = 2 



I677 

2(<7j_ + <7| 

<7j_ +2(7|| 
1+2(711 



(l+2(7 ± + 6(7||) 2 
8(1+2(7,1) 



(53) 
(54) 

(55) 

(56) 



In the isotropic case, we get the growing kinematic solution: 
7 = (3 /4)7 and 5 = 3/2 (see §3.1). This growth can be quenched 
by making velocity anisotropic so that it varies more slowly 
along the folds than across: the values of a± and (7|| then drop 
compared to the isotropic case, and so does the growth rate 7 
until the dynamo is shut down. Thus, saturation can be achieved 
by anisotropizing the statistics of the velocity gradients. 

The parallel gradients of the velocity field that occur in 
the numerators a± and <7|| [eqs. (49-50)] represent stretch- 
ing motions, while the perpendicular gradients that contribute 
to 7^ [eq. (48)] correspond to the two-dimensional mixing 
(interchange-like) motions that shuffle and bring the direction- 
reversing magnetic fields sufficiently close together for them to 
annihilate resistively. The decrease in the values of a± and a\\ 
means that the comparative strengths of mixing and stretching 
change in favor of mixing. The steady solution is thus a result 
of a balance between weakened stretching and two-dimensional 
mixing of the folded fields by the partially two-dimensionalized 
random flow. 

Note that the mixing rate 7^ may or may not be changed by 
the nonlinearity, with dynamo saturation achievable either way, 
provided that the stretching rates weaken sufficiently. What 
happens to 7^ is important because it enters the expression for 
k v [eq. (54)] and thereby determines the shape of the saturated 
spectrum. It can be seen in § 5.2 that whether mixing is sup- 
pressed is a key issue for the large-Re dynamo. For now, we 
assume that j± is unchanged and see if we can get a reasonable 
model of saturation for our viscosity-dominated runs. 



4.3.2. Solution of the Model 

The condition for saturation is 7 = [see eq. (55)]. This, 
however, does not uniquely fix cr± and cry , but only relates them 
via a quadratic equation. Thus, to get a specific solution, we 
must make a further assumption about cr± and (7|| . One plau- 
sible argument is that the anisotropization of the velocity gra- 
dients occurs in such a way that all modes with k above some 
self-consistently chosen cutoff k s (the stretching wavenumber) 
are completely two-dimensional (have = 0), while the mo- 
tions with k < k s are still isotropic (this is certainly a sensible as- 
sumption in the case of large Re, see Schekochihin et al. 2004b, 
and § 5.2). Then the reduction in the integrals in eqs. (49-50) is 
due to the shrinking of the integration domain in k space, but the 
integrands are still isotropic [k (i) = K (l) (fe), K (a) = for k < k s ], 
which implies that (7j_/(7|| = 4 as in the isotropic case. 

With this assumption, the dynamo saturates (7 = 0) when 
a± = 4(7|| ~ 0.078. The spectral exponent [eq. (56)] is s ~ 0.23, 
which, however, is clearly a model-dependent number. 

4.3.3. Matching with Simulation Results 

The solution (53) is valid in the limit k v ^>k*~ k v (Pr m » 1), 
but, as we are about to see, convergence in Pr m is only loga- 
rithmic. Numerical solution of eq. (47) shows that scale sep- 
arations of several decades are required for the scaling to be 
discernible. No such scale separations are achievable in direct 
numerical simulations. Therefore, in order to make a meaning- 
ful comparison, we must consider model-predicted spectra for 
the case of finite separations between the resistive cutoff k v and 
the infrared cutoff k*, where the boundary condition is imposed. 
Seeking steady solutions of eq. (47), we get 



M(k) = (const) k s K i0 . 



'(C_L,0-||) 



(k/kjj), 



(57) 



1 /2 

where D(a±,ai\) = [8(7/7^/(1 + 2(7||)] and k v , 7, s are de- 
fined in eqs. (54-56). The equation for a± and <7|| is obtained 
from the zero-flux boundary condition: the expression in square 
brackets in eq. (47) must vanish at k = k* (i.e., no energy passes 
through the left boundary in k space). This gives 

k* (k„ \ l+2(7 ± + 6(7|| (K\ 

Together with the assumption <r± = 4<7|| , this is a transcendental 
equation for <j\\ , which can be solved numerically for any fixed 
set of parameters k*, rj, -f±. 

Note that in the limit of k* <C k v , eq. (58) reduces to 
sin[i>ln(A: >t /2A:^)] = 0, whence 



7 = 



d+2(7||)^ 2 

8[ln(*./2ifc,)] 



-7-L 



O 



1 



[lnPr, 



1/2-1 



(59) 



As Pr m — > 00, we have 7^0, converging logarithmically. This 
is the asymptotic solution discussed in §4.3.1 and §4.3.2. 

To make a comparison with our numerical results, we set k* = 
ko = 2tt, the box wavenumber, and use the following (isotropic) 
expression for the mixing rate: 7^ = (6/5)7 = c ( e /^) 1,/2 > where 
v is viscosity, e = 1 is the mean forcing power [see eq. (4)], and c 
is an adjustable parameter. We have compared the model solu- 
tions with the (normalized) spectra obtained in our viscosity- 
dominated runs. The adjustable parameter c was the same 
(c = 0.3) for all cases. The entire sequence is well fitted by 
the model (except at k/2-K = 1,2, where boundary conditions 
and discretization effects are important): see Fig. 20. 
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FIG. 20. — Magnetic-energy spectra (normalized by the total magnetic en- 
ergy) in the saturated state. The thin red lines are spectra predicted by our 
model. For run S 1 , the model does not give energy growth, but Pr m in this case 
is too small for any asymptotic consideration to be applicable. 

Considering the simplicity of our model, its ability to pre- 
dict nonasymptotic numerical spectra so well is remarkable. 
It shows that numerical results are consistent with magnetic- 
energy spectra that are peaked at the resistive scale and would 
have a very flat positive spectral exponent in the interval k v -C 
k -c k n if sufficiently large scale separations were resolved. 

4.4. Intennittency 

In the saturated state, intermittency is reduced, implying 
tighter packing of the system domain by the magnetic field 
(Fig. 11a). It is not a surprising result: nonlinear back- 
reaction imposes an upper bound on the field growth, and, once 
the strongest fields in the dominant substructure saturate, the 
weaker ones elsewhere have an opportunity to catch up. The pdf 
of the saturated field is exponential at large B (Fig. 13, cf. Bran- 
denburg et al. 1996; Cattaneo 1999). In the saturated state, the 
kurtosis does not depend on Pr m (Fig. 1 lb). Furthermore, while 
B rms is not Pr m -independent for finite values of Pr m (Fig. 14), 
the pdf of B/B lms is the same for all Pr m . 

For B < B rms , the pdf is power-like, P(B) ~ B 1 - 7 - 1 - 8 , and un- 
changed from the kinematic regime (see § 3.4.3 and Fig. 12b). 
This tail is made up of the weak fields in the bends of the folds. 
As we argued in §4.1 in the context of the curvature distri- 
bution, these fields remain quasikinematic, with no significant 
change in their statistics compared to the kinematic regime. 

4.5. Summary of the Saturated State 

We now summarize what we have learned about the saturated 
state of the dynamo in the viscosity-dominated regime: 

1. The folded field structure inherited from the kinematic 
regime is preserved: magnetic fields remain organized 
in flux sheets (or ribbons) with alternating field direction 
and field strength and field-line curvature are anticorre- 
lated. The scale separation between the fold length (flow 



scale) and the direction-reversal scale (resistive scale) is 

pj.l/2 

2. The velocity field in the subviscous scale range (k > k„) 
develops a k~ A spectrum, which is induced by the subvis- 
cous fluctuations of the Lorentz tension force and plays 
no part in the nonlinear back-reaction. 

3. The saturated magnetic-energy spectrum is flatter than 
in the kinematic case, but the bulk of the energy is still 
at the resistive scale. Accessible numerical resolutions 
do not allow us to measure the spectral exponent. How- 
ever, the numerically obtained spectra are in an excellent 
agreement with those predicted for the same values of 
Pr m by a simple saturation model, which, in the limit of 
very large Pr m , would give a very flat positive exponent. 
This model is based on the idea that the nonlinear sat- 
uration is achieved via anisotropization of the velocity- 
gradient statistics with respect to the local direction of 
the folds. The anisotropization consists in a reduction 
of the stretching component of the flow compared to its 
mixing component. The saturated spectra result from a 
balance between the weakened stretching and the mixing 
of the fields by quasi-2D interchange-like motions. 

4. The fields in the saturated state are somewhat more 
tightly packed than in the kinematic regime, but are still 
intermittent: instead of a lognormal pdf, they have an ex- 
ponential one. Unlike in the kinematic case, the volume- 
filling fraction is independent of Pr m . 

5. SMALL-SCALE DYNAMO IN TURBULENCE WITH LARGE 
REYNOLDS NUMBERS 

Up to this point, we have not, in fact, attempted to study 
the turbulent dynamo in the proper sense. Because of the large 
viscosity, our velocity was spatially smooth and single-scale. 
The only feature that distinguished our model from synthetic 
dynamos in prescribed flows (such as those reviewed in Chil- 
dress & Gilbert 1995) was that our velocity was random in 
time (because of the random forcing), homogeneous in space 
(had no fixed spatial features such as stagnation points), and 
subject to the physical nonlinear back-reaction associated with 
the Lorentz-tension term in eq. (1). The main motivation for 
dwelling on this regime at such length was that a parameter 
scan in Pr m could be afforded within accessible resolutions. 
However, if the small-scale dynamo theory is to make state- 
ments about physical reality, its ultimate goal must be to under- 
stand small-scale dynamo in turbulent velocity fields with large 
Reynolds numbers. In other words, we must attack the problem 
of isotropic MHD turbulence. 

This problem is essentially distinct from the more studied 
anisotropic case in which a strong mean field is imposed on the 
system. This latter case is somewhat better understood, mostly 
because we have what probably is a correct heuristic view of 
it as a turbulence of interacting Alfven-wave packets (due to 
Kraichnan 1965), just as our intuition about Kolmogorov tur- 
bulence is based on the idea of the constitutive fluid motions 
as "eddies" [note, however, that the more recent and increas- 
ingly widely accepted phenomenology of Goldreich & Sridhar 
(1995) with its idea of critical balance has blurred the distinc- 
tion between Alfven-wave packets and eddies]. 

It is the presence of the mean field that makes Alfven waves 
possible. If no mean field is imposed, it is tempting to assume 
that we are still dealing with Alfven waves, except now waves 
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FIG. 21. — Growth and saturation of the magnetic energy for (a) Pr m = 1 runs (run al saturates around t = 50, not shown) and (b) run B. 



at small scales see the outer-scale magnetic fluctuations as the 
local mean field, along which these waves can propagate. How- 
ever, in order for such a scheme to work, the outer-scale fluctua- 
tions must be more energetic than the small-scale ones, i.e., the 
bulk of the magnetic energy must be self-consistently kept at 
large (box) scales by some nonlinear mechanism. We are about 
to see that this does not happen (§ 5.1.2). 

Resolution constraints do not permit a meaningful parameter 
scan in the large-Re and large-Pr m regime, so we limit our dis- 
cussion to two representative runs at the 256 3 resolution: run A 
is a traditional Pr m = 1 run with v = 5 x 10" 4 , while run B has 
Pr m =10 and u = 2x 10" 3 . Let us define the Reynolds number 
as 



Re = 



(60) 



where ko = 2n is the box wavenumber. Then in the saturated 
(kinematic) stage we have Re ~ 390 (450) for run A and Re ~ 
80(100) for run B. The corresponding values of the Taylor- 
microscale Reynolds number defined by 



«« 2 >/3) 1/2 A 



Rex = 



1/2 



<« 2 > 



1/2 



vk\ 



(61) 



(see Frisch 1995, and eq. (16)) are Re A ~ 155(116) for run A 
and Re A — 45 (50) for run B. In a few instances, we shall com- 
pare run A with two other Pr m = 1 runs (al and a2; see Table 1) 
that have lower Re. None of these runs are truly asymptotic in 
either Pr m or Re, so we do not attempt to extract scalings and 
all our conclusions are tentative. 

Note that even in the case of Pr m = 1 (Fig. 21), the resolu- 
tion constraints for isotropic MHD are more severe than for 
purely hydrodynamic simulations because (1) the resistive scale 
is smaller than the viscous scale (see Fig. 22a and Fig. 23a) and 
(2) in order to achieve statistical steady state and collect ade- 
quate statistical information, simulations must be run for tens 
of box-crossing times, not just one or two as in the hydrody- 
namic case. The highest-resolution simulation of incompress- 
ible isotropic MHD to date is a 1024 3 one with Re ~ 960 and 



Pr m = 1 by Haugen et al. (2003), who were able to run it for a 
few box-crossing times. It appears, however, that even larger 
resolutions and Pr m ^> 1 runs are needed to make statements 
about the truly asymptotic case. 

In §5.1, we describe the growth, saturation, and spectra of 
the magnetic energy in runs A and B on a purely factual level. 
A theoretical discussion is given in §5.2. In §5.3, we return 
to the numerical results to test our theoretical arguments and 
to make some tentative conclusions (§5.3.5). Intermittency is 
studied in § 5.4. The summary of the main results is in § 7. 

5.1. Growth, Saturation, and Energy Spectra: Facts 
5.1.1. The Kinematic Regime 

The small-scale kinematic dynamo in the large-Pr m limit de- 
pends very little on the particulars of the velocity field (as long 
as this velocity is three-dimensional and sufficiently chaotic). 
As seen below, all indications are that our runs A and B are 
nonasymptotic cases of the large-Pr m regime; i.e., the same 
stretching and folding mechanism is responsible for the small- 
scale-field generation. In their kinematic stage, these runs are, 
therefore, very similar to the viscosity-dominated simulations. 

Both runs exhibit exponential magnetic-energy growth in 
the kinematic regime (Fig. 21). In Kolmogorov turbulence, 
small eddies turn over faster than the large ones, so dominant 
stretching is done by the viscous eddies. We expect, therefore, 
that the growth rate is 7 <~ r rms when Re » 1, where r rms = 

1 /2 

((|Vu| 2 )/15) is the rmsrateof strain. Since j/(|Vu| 2 ) = e= 1 

[eq. (4)], we have r rms ~ Re 1/2 r^ ? , where r box = (u 2 ) 1/2 is the 
box crossing time (the box size is 1). In principle, to check 
these statements, we need a Re parameter scan at constant large 
Pr m . While this is not possible at current resolutions, Pr m = 1 
runs can be used to establish qualitatively that 7 increases with 
Re. In Table 2 we give values of 7 (least-squares fit; see also 
Fig. 21a), T ms , and r^J x for run A and for two extra runs with 
Pr m = 1 and smaller Re. For comparison, we also include the 
data for runs B and S4. Even though the numerical values of 7 
are substantially smaller than r rms (and might appear, mislead- 
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FIG. 22. — Energy spectra for (a) run A and (b) run B in the kinematic regime and in the saturated state. The magnetic-energy spectra in the kinematic regime are 
normalized by the total magnetic energy, then averaged and rescaled by the total magnetic energy of the saturated state (for comparison with the saturated spectra). 



ingly, to be closer to T bc J x ), it is clear that both r rms and 7 grow 
with Re, while r^J x stays the same. The smallness of the values 
of 7 relative to the rms rate of strain should probably be at- 
tributed to the resistive correction being order-unity for Pr m = 1 
and moderate values of Rm in our runs [see eq. (12)]. 

It is fair to observe that the increase of 7 with Re that we 
see has, in fact, two causes: the growth of the rate of strain 
with Re and the increasing value of Rm. Indeed, in run al, 
Rm ~ 100 is quite close to the critical value Rm c , below which 
there is no dynamo: this is demonstrated by run aO (Table 1), 
which has Rm ~ 52 and decays [this implies Rm c e (52, 100), 
which is the same estimate as that for the viscosity-dominated 
runs; see § 3.1]. The effect of finite Rm on the value of 7 will 
saturate at large Rm, while the large-Re asymptotic should be 
7 <~ Re 1 / 2 . The values of Re in our runs are not large enough 
for us to be able to make statements about scalings. Simulations 
at larger Re appear to be consistent with the scaling 7 ~ Re 1//2 
(N. E. Haugen 2003, private communication). 

In both runs A and B, the magnetic energy, initially at the 
box scale, is quickly piled up the resistive scale, just as in the 
viscosity-dominated case. The evolution of the spectra is again 
self-similar (analogous to Fig. 4a). The normalized spectra av- 
eraged over the kinematic regime are shown in Figure 22. 

5. 1 .2. The Saturated State 

The dynamo eventually saturates and a statistical steady state 
is achieved (Fig. 21). The kinetic energy (u 2 ) in this regime is 
smaller than in the purely hydrodynamic case (B = 0). The mag- 
netic energy (B 2 ) is smaller than (u 2 ) by a factor of 6 in run A 
and by a factor of <~ 2.3 in run B (Table 1). The saturated spec- 
tra (Fig. 22) exhibit no scale-by-scale equipartition between the 
kinetic and magnetic energies even at Pr m = 1 . 

In both runs A and B, the kinetic-energy spectrum in the 
purely hydrodynamic regime has a very short segment consis- 
tent with the Kolmogorov A:" 5 / 3 scaling (this is somewhat better 
pronounced for run A). In the nonlinear saturated MHD state, 
the spectrum is steeper: in each of the runs, there is a segment 



best fitted by A;" 7 / 3 . It is not clear if this scaling (if it is re- 
ally there) survives in the large-Re limit or is simply a limited- 
resolution effect (spectra steeper than Kolmogorov sometimes 
appear in hydrodynamic simulations with too large viscosity; 
see theoretical discussion in Frisch 1995). If this scaling is, 
in fact, asymptotic in Re, it is a feature of the fully-developed 
isotropic MHD turbulence for which, at present, we have no 
theoretical explanation (see some further observations in § 5.2). 
What is certainly true is that (1) the kinetic energy is dominated 
by the large (box) scales and (2) the characteristic "eddy" time 

1 /2 

<~ [A: 3 £■(&)] peaks at the viscous scale. 

The magnetic-energy spectrum is shifted toward small (resis- 
tive) scales, with a substantial excess of magnetic over kinetic 
energy there. There is no discernible scaling. Note that the 
shapes of the magnetic-energy spectra in runs A, B and S1-S6 
(Fig. 20) are not very different. Thus, it appears that there is no 
inverse cascade of magnetic energy and that the characteristic 
scale of the magnetic field is the resistive scale. 

The fact that numerical simulations do not support equiparti- 
tion was first explicitly emphasized by Maron et al. (2004), even 
though many prior numerical studies had found similar spectra 
(Meneguzzi et al. 1 98 1 ; Kida et al. 1 99 1 ; Miller et al. 1 996; Cho 
& Vishniac 2000; Chou 2001). However, there was a reluctance 
to accept that numerical evidence was at odds with the equipar- 
tition hypothesis and, therefore, with the standard Alfvenic phe- 
nomenology. In order to extract scalings, a popular device (due 
originally to Kida et al. 1991) has been to consider the spectrum 
of the total energy, E(k) + M(k). In our opinion, this choice of 
diagnostic obscures rather than clarifies the situation. Indeed, it 
only makes sense to add the spectra together if there is a sym- 
metry between u and B. For the spectra that obtain for Pr m = 1 
and at limited resolutions, adding E(k) and M(k) gives one what 
might be a false impression that a Kolmogorov scaling actually 
emerges, as shown in Figure 22a. However, without equipar- 
tition, there is no obvious reason why asymmetric kinetic and 
magnetic energy spectra should add up to something that has a 
Kolmogorov scaling. 
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One could argue that all current simulations are nonasymp- 
totic and that, were the resolution to be increased dramatically, 
equipartition and Kolmogorov scaling would emerge at hith- 
erto inaccessible small scales. While this cannot be ruled out, 
the indications so far are not encouraging. Haugen et al. (2003, 
2004) performed simulations at resolutions ranging from 64 3 
to 1024 3 (with Pr m = 1 and maximum Re = 960). In none of 
their runs were the kinetic and magnetic energies in scale-by- 
scale equipartition, with magnetic energy always shifted to- 
ward the small scales. Furthermore, they found that the Kol- 
mogorov scaling that seemed to hold at 256 3 (cf. Fig. 22a) ac- 
tually got worse at 1024 3 , with magnetic energy bulging up at 
small scales. They attributed this problem to some form of the 
bottleneck effect, i.e., essentially, to nonlocal interactions in k 
space. Small-scale dynamo is, in fact, just that: a nonlocal in- 
teraction between velocity and magnetic fields. 

A Kolmogorov scaling for the total energy was also reported 
by Biskamp & Miiller (2000) who simulated decaying isotropic 
MHD turbulence with Pr m = 1 at resolutions of 512 3 . It is, how- 
ever, unlikely that decaying-MHD results apply to the forced 
case. Indeed, the kinetic and magnetic energy spectra that 
Biskamp & Miiller (2000) found in their self-similarly decay- 
ing regime were quite different from those for the forced case: 
while again there was no equipartition, it was the magnetic en- 
ergy that dominated the kinetic one at large scales (we presume 
that this is due to large-scale MHD equilibrium structures). 

Thus, on the balance of the available numerical results, while 
we cannot definitively claim that the equipartition theory is 
wrong, it is certainly true that there is no evidence in support 
of it and a fairly strong case against. 

5.2. The Physics of Saturation: Selective Decay, Alfven 
Waves, and Mixing 

Let us consider what are the reasonable physical scenarios of 
nonlinear transition from field growth (dynamo) to saturation 
and what they imply about the final state. 

5 .2. 1 . The Intermediate Nonlinear Regime 

The folded field structure implies the criterion for the onset 
of nonlinearity: the back-reaction is controlled by the Lorentz 
tension force B • VB <~ k\\B 2 , which only depends on the parallel 
gradient of the field and does not know about direction rever- 
sals (kii <~ k u ). Balancing B • VB <~ u • Vu, we find that back- 
reaction is important when magnetic energy becomes compara- 
ble to the energy of the viscous-scale eddies. Eddies at scales 
above viscous are still more energetic than the magnetic field 
and continue to stretch it at their (slower) turnover rate. When 
the field energy reaches the energy of these eddies, they are 
also suppressed and the dominant remaining stretching is done 
by yet larger and slower eddies. Consider a time t when the to- 
tal magnetic energy (B 2 )(t) is less than its final saturation value 
but larger than the energy of the viscous-scale eddies. Let the 
stretching scale l s (t) be the scale of the eddies whose energy is 
equal to (B 2 ), i.e., 



(B 2 ). 



(62) 



All eddies smaller than l s have had their stretching component 
suppressed, while all the larger eddies are slower, so, at time t , 
the effective rate at which the field is being stretched is uijl s . 
Using this and eq. (62), we estimate 



ui. 



dt 



(B ) ~ ^ (B } ~-±~e = const, 



(63) 



where e is the Kolmogorov energy flux. Equation (63) implies 
(B 2 )(t) ~ et, i.e., there is an intermediate period of nonlinear 
evolution between the onset of nonlinearity and saturation, dur- 
ing which the magnetic energy grows proportionally to time 
rather than exponentially. The force balance underpinning the 
back-reaction by the folded fields on the eddies of size / s is 



u • Vu <~ — - 

h 



'B-VB-/t|,(B 2 ), 



(64) 



so that <~ l~ l , i.e., the folds are elongated to the stretching 
scale (for a slightly more quantitative argument in favor of fold 
elongation, see Schekochihin et al. 2002d). 

The key question is by how much l s can exceed the viscous 
scale Iv, for it is the saturated value of l s that will determine 
the saturated magnetic energy and its spectral distribution. Sat- 
uration is achieved when the nonlinear growth according to 
eq. (63) is balanced by diffusion: 

T 72 ■ (65) 

Therefore, in order to determine 4, we need to know what hap- 
pens to the resistive scale l v in the nonlinear regime. There are 
two main possibilities (subject to the caveats in § 6.1). 

5.2.2. Saturation with Efficient Mixing 

If the velocities at I < l s act as interchange-like motions in 
the plane transverse to the local direction of the folds and, con- 
sequently, mix the field lines in a quasi-2D fashion, then the 
resistive scale is determined by the characteristic timescale of 
these motions at the viscous scale, 



I 2 ~ / ' 



(66) 



Comparing this with eq. (65), we see that l s cannot grow larger 
than l v by more than a factor of order unity. In this case, the 
motions at or just above the viscous scale are anisotropized, so 
mixing is enhanced and wins the competition with stretching, 
leading to saturation. The magnetic energy saturates at a value 
determined by the energy of the viscous-scale eddies. Thus, 
the signature of this type of saturation would be the following 
dependence of the saturated magnetic energy on Re: 

o2\ 



B z 



Re 



■ l/2 (u 2 ) 



(67) 



5.2.3. Saturation with Inefficient Mixing 



The second possibility is that mixing by motions at I < l s 
is inefficient, so that both stretching and mixing of the field 
is mostly done by the stretching-scale eddies. It is then the 
turnover time of these eddies that determines the resistive scale 
of the magnetic fields, giving [from eq. (65), (62), and (63)] 



In 



1/2 



(Jit) 



1/2 



(68) 



Physically, this expression describes selective decay of the 
magnetic energy at scales too small to be sustainable by the 
weakened stretching (or mixing). In this scenario, the magnetic 
energy is free to grow to the equipartition level, 

(B 2 )~(u 2 ), (69) 

at which point l s <~ Iq and the resistive scale is determined by 
the stretching rate of the outer-scale eddies: 



L~(^) -Rm-^o-Re^Pr- 1 /^, 



1/2 



(70) 
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Evolution of characteristic wavenumbers (denned in § 3.2.1) for (a) run A and (b) run B. Note that for Kolmogorov turbulence, k\ ~ Re l l A k v . 



FIG. 23.— 

i.e., it is larger by a factor of Re 1 / 4 than in the kinematic regime. 
Note that the condition for l v to be distinquishably smaller than 
l v in a numerical experiment is, therefore, 

Pr m >Re 1/2 > 1, (71) 

which is impossible to satisfy at present resolutions. Since the 
folds elongate so that £11 <~ l~ l <~ Iq 1 (§5.2.1), the separation be- 
tween the direction-reversal scale and the fold length increases: 

^-^Re^Pri/ 2 , (72) 

K|| l n 

which is a factor of Re 1 7/2 larger than in the kinematic regime. 

If the velocities at I < l s are not interchange-like mixing mo- 
tions, what are they? There is another type of motion that does 
not amplify the magnetic field and can, theoretically, be present 
below the stretching scale: Alfven waves can propagate along 
the folded direction-reversing magnetic fields (Fig. 24). These 
waves are damped by viscosity and make a contribution to both 
the kinetic and magnetic energy spectra. Their dispersion rela- 
tion is lu 2 = kk : bb(B 2 ) for k v <C k <C k\\ ~ l~ l (Schekochihin 
et al. 2002d). This interval is non-empty provided that l s 3> /„. 
Therefore, mixing must be inefficient in a turbulence made up 
of such waves; otherwise, / s cannot grow sufficiently far above 
l v for them to exist in the first place. 




FIG. 24. — Alfven waves along folded fields: a sketch. 

Detecting these Alfven waves in numerical simulations is a 
considerable challenge both because of large scale separations 
that must be resolved and because designing the appropriate di- 
agnostics is not straightforward. The simulations we present 



here are not appropriate for such a task. However, some of the 
physical possibilities we have anticipated can be tested. 

5.3. Numerical Evidence: Facts and Interpretations 

5.3.1. Saturated Energy Levels 

The key challenge now is the resolution of the dichotomy de- 
scribed in § 5.2.2 and § 5.2.3: saturation with magnetic energy 
comparable to the viscous-eddy energy [eq. (67)] or to the total 
energy of the turbulence [eq. (69)]. The most direct approach 
would be to look at the scaling of the ratio (B 2 ) / (u 2 ) with Re at 
constant large Pr m , but this is not yet possible at current resolu- 
tions. At Pr m = 1, (B 2 ) /(u 2 ) increases with Re (Table 2). In the 
higher-resolution Re scan with Pr m = 1 produced by Haugen et 
al. (2003), (B 2 )/(u 2 ) appeared to saturate at larger Re. Similar 
results were reported in earlier numerical studies (Kida et al. 
1991; Kleva & Drake 1995; Miller et al. 1996; Archontis et al. 
2003b). Thus, at least for Pr m = 1, numerical results do not sup- 
port the scaling (67), suggesting that it is the second possibility 
(§ 5.2.3) that is realized. 

5.3.2. The Intermediate Nonlinear Regime 

We now examine the evolution of our runs A and B from the 
kinematic to nonlinear regimes to see if we can identify any of 
the features of the saturation scenarios proposed in § 5.2. 

We notice immediately that, between the kinematic stage of 
exponential growth and the saturated statistically steady state, 
an intermediate nonlinear stage during which the magnetic en- 
ergy grows slower than exponentially is manifest in run A and, 
to a lesser extent, in run B (Fig. 21). This intermediate regime 
was not present in the viscosity-dominated runs and is clearly a 
large-Re effect in agreement with the argument in § 5.2.1. 

5.3.3. The Folded Structure 

Employing the same diagnostics as in § 3 and § 4, we find 
that the folded structure appears to survive in runs A and B de- 
spite the moderate values of Pr m . 

The time histories of the characteristic wavenumbers (de- 
fined in §3.2.1) are shown in Figure 23. In the kinematic 
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Fig. 25. — Curvature pdf's for (a) run A and (b) run B during the kinematic regime and in the saturated state. The slopes corresponding to K 1 and K n l 1 [see 
eq. (32)] are given for reference. 



regime, k n 



^Bxj 3> k\\- During the intermediate nonlin- 



ear regime, h\ decreases further and becomes comparable to 
the box wavenumber. The resistive wavenumber k lms also de- 
creases, but by a lesser amount, which is consistent with the 
ideas of selective decay, fold elongation and increased scale 
separation in the folds in the saturated state set forth in § 5.2.1 
and § 5.2.3 [see eq. (72)]. The fold elongation is a new effect: 
it could not be seen in the viscosity-dominated runs (§4), for 
which there was no separation between the viscous and box 
scales. 

All other features associated with the folded structure that 
we reported in § 3.2 and § 4.1 for the viscosity-dominated runs 
carry over to runs A and B: the magnetic-field strength and cur- 
vature are anticorrelated (e.g., the correlation coefficient Ck,b, 
eq. (26), is within 6% of -1, see Table 1) and the pdf's of cur- 
vature (Fig. 25) are again concentrated at velocity scales. The 
power tail appears to be somewhat steeper than K~ n l n predicted 
for the large-Pr m case (the steepening is more pronounced for 
run A). Interestingly, the scaling Pk(K) <~ K at K <C k v (to the 
left of the peak) suggested by eq. (32) seems to be in evidence 
for the kinematic pdf's. The elongation of the folds in the sat- 
urated state compared to the kinematic regime is manifested by 
the peak of the curvature pdf shifting to the left. 

5.3.4. The Tension Spectra 

Let us now look at the magnetic back-reaction B • VB as a 
function of k: the tension spectra for runs A and B are plot- 
ted in Figure 26. For run B, the tension spectrum is flat up to 
the outer (box) scale; for run A, it has a plateau that starts at 
k/2ir <~ 3, which should probably also be interpreted as corre- 
sponding to the outer scale. Based on the argument of § 3.2.2, 
such flat tension spectra are consistent with magnetic fields that 
are folded with the characteristic fold length at the outer scale. 
In a further confirmation of the consistency of this view, the 
relation T(k) ~ kjM 4 (k) [eq. (23)] roughly holds (see Fig. 26). 

Thus, the tension spectra in runs A and B have similar prop- 
erties to those in the viscosity-dominated case (§4.2) except 



the tension force B • VB is no longer passively compensated 
by the velocities as in eq. (43): Fig. 26 shows that the relation 
(45) is only satisfied at very small scales. Instead, in the satu- 
rated state, we expect the nonlinear balance (65) to hold at the 
stretching scale l s with velocities at / > / s remaining essentially 
hydrodynamic and velocities at I < l s modified so that they can- 
not significantly affect the folded structure, i.e., either mix or 
stretch the magnetic field. Since eq. (65) is only an order-of- 
magnitude statement, it is not meaningful to compare the actual 
values of u • Vu and B • VB. However, it is still instructive to 
compare the behavior of T(k) with that of k 3 E(k) 2 , a quantity 
that roughly measures the spectrum of the hydrodynamic inter- 
action term u • Vu. In Figure 26, this quantity is plotted in both 
the purely hydrodynamic case and the saturated MHD state 
[note that the absolute magnitude of k 3 E(k) 2 has been rescaled 
in Figure 26 and should not be compared to that of T(k)]. In- 
terestingly, the intervals of flat T(k) for both runs coincide with 
the intervals where there is a substantial suppression of k 3 E(k) 2 
in the MHD case compared to the hydrodynamic case. These 
intervals correspond to the steeper-than-Kolmogorov scaling of 
the kinetic-energy spectra noted in § 5.1.2. While we reiterate 
the caveat made in § 5.1.2 that this scaling might be an artifact 
of limited resolution, it is just possible that we are dealing with 
a genuine spectral signature of the nonlinearly modified fluid 
motions in the interval k$ < k < k v . 

Based on these plots, we hypothetically identify the stretch- 
ing wavenumbers for our runs: for run A, k s /2n ~ 3; for run B, 
k s /2TT ~ 1. The resistive and viscous cutoffs for these runs are 
all around k/2n ~ 10 and hard to tell apart in the saturated state, 
which is consistent with the estimates in §5.2.3 and Pr m > 1 
[note that even run B does not really satisfy the condition (71) 
necessary for k v and k v to be distinguishable]. 

5.3.5. Tentative Conclusions 

At the fairly low values of Re that our runs have, we can- 
not test scalings such as eq. (67), (69), (70), or (72) so as to 
make a definite distinction between the two scenarios of satura- 
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tion (§ 5.2.2 vs. § 5.2.3). On the qualitative level, however, our 
results appear to favor the second scenario: 

• increase of the ratio (B 2 ) /(u 2 ) with Re, 

• the presence of an intermediate nonlinear regime, 

• the evolution of the characteristic wavenumbers (selec- 
tive decay and fold elongation), 

• the evidence of the folded structure up to the box scale 
(curvature pdf's, anticorrelation between curvature and 
field strength, flat tension spectra), 

• the nonlinear modification of the kinetic-energy spectra 
up to the box scale 

are all consistent with the conclusion that the stretching scale 
in saturation is comparable to the outer scale, not the viscous 
scale, 13 as suggested in § 5.2.3. In this second scenario of sat- 
uration, the velocities at scales below Z s ~ Iq, do not destroy 
the folded direction-reversing magnetic fields. The magnetic 
fluctuations are of two kinds: the folds of length Iq with direc- 
tion reversals at l v [eq. (70)] and (hypothetically) Alfven waves 
propagating along the folds with k in the inertial range (§ 5.2.3). 
The dynamo is saturated at an energy that is comparable to (w 2 ). 

The mechanism of this saturation should be entirely con- 
trolled by what happens at the outer scale. Since the motions 
at this scale are directly excited by the forcing, the nonlinear 
dynamics there may not be totally dissimilar from the nonlin- 
ear dynamics in the viscosity-dominated case. In other words, 
we might conjecture that saturation results from a competition 
of weakened stretching and enhanced mixing at the outer (forc- 
ing) scale (the hypothesis that the dynamo saturation is con- 
trolled by the outer timescale was first explicitly formulated by 
Maron et al. 2004). 



In §4.3, we constructed a quantitatively successful model 
of saturation in the viscosity-dominated case by assuming that 
velocity gradients are partially anisotropized leading to weak- 
ened stretching but unsuppressed mixing. Since the viscosity- 
dominated runs have Re <~ 1, the stretching scale is the box 
scale, / s ~ Iq, which is also the viscous scale, Iq ~ l v , so the 
success of this model cannot be construed to favor either of the 
saturation scenarios outlined in § 5.2. The fact that the mixing 
rate did not have to be suppressed in order for correct results 
to be obtained is consistent with both the scenario of efficient 
mixing of the field by the eddies at the viscous scale l v and the 
idea that the saturated state is controlled by the balance between 
mixing and stretching at the box scale Iq. 

The reader is referred to § 6.1.2 for some further caveats. 

5.4. Intermittency 

We complete our standard set of diagnostics by looking at 
the pdf's of the magnetic-field strength in runs A and B. 

In the kinematic regime, their evolution is again self-similar 
with the pdf's of B/B lms collapsing onto stationary profiles 
(similarly to Fig. 12). These profiles (averaged over the kine- 
matic stages of the runs) are shown in Figure 27. The lognormal 
fit [eq. (42)] works even better for these runs than it did for the 
viscosity-dominated case in §3.4.3 (Fig. 13). 

The intermittency is reduced in the saturated state: the evo- 
lution of the normalized moments is analogous to the viscosity- 
dominated case (Fig. 11a) with the values of kurtosis in the sat- 
urated state still super-Gaussian, but substantially smaller than 
in the kinematic regime (see Table 1). On a semilog plot, the 
pdf's in the saturated state of both runs (Fig. 27) have straight 
segments in the region of strong fields (exponential tails). How- 
ever, at even larger B, there appears to be a knee followed by 
steeper decay. To a lesser extent, the same feature was also 
present in the viscosity-dominated case (Fig. 13). We believe 



13 In principle, it is not excluded that the stretching scale l s saturates at some Re-dependent value intermediate between Iq and l v . In fact, one might consider the fact 
that k s appears to be somewhat larger than the box wavenumber in run A suggestive of such an outcome. Much higher resolutions are necessary to test this possibility. 
However, the scale invariance of the inertial range makes saturation with Iq -C Z s -C l v a priori unlikely. 
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FIG. 27. — The pdf's of B/B r ms for (a) run A and (b) ran B averaged over the kinematic stages of these runs and the corresponding lognormal profiles [eq. (42)] 
with the same D. The pdf's in the saturated state are also given. 



that this feature might be due to resolution effects and/or insuf- 
ficient averaging (as these values of B correspond to very rare 
events), so the physical pdf at large B is exponential. 

Note that exponential pdf's of the magnetic-field strength 
were reported previously in Pr m > 1 simulations of compress- 
ible convection with rotation (Brandenburg et al. 1996) and of 
Boussinesq convection (Cattaneo 1999). Thus, exponential tails 
of the field-strength pdf appear to be a fairly universal feature in 
MHD turbulence (at least away from the boundaries; see Catta- 
neo 1999; Thelen & Cattaneo 2000). 

6. DISCUSSION 
6. 1 . Are All Small-Scale Dynamos Large-Pr m Dynamos? 

The fundamental mechanism of magnetic-energy genera- 
tion by three-dimensional chaotic and/or random flows (due to 
Zel'dovich et al. 1984) is based on the idea that the structure of a 
locally uniform rate-of-strain tensor in three dimensions allows 
for a winning configuration of magnetic fields: the field aligns 
with the stretching direction of the flow but varies (reverses) 
along the "null" direction, so that compression cannot cancel 
the effect of stretching (see § 3.3). Thus, the folded structure, 
which has been the leitmotif of this paper, is an essential ingre- 
dient in our understanding of the very existence of the dynamo. 

The Zel'dovich et al. mechanism is universal in the sense that 
the particulars of the flow matter little. However, two proper- 
ties are essential: (1) the flow must be spatially smooth and (2) 
the magnetic field must be able to vary at scales much smaller 
than the scale of the flow. Only if these conditions are satis- 
fied, can one think of the flow as locally linear and make the 
above inferences about the growing field's structure. The tur- 
bulent dynamo with large Pr m satisfies these conditions and so 
do all deterministic chaotic dynamos (see reviews by Childress 
& Gilbert 1995; Ott 1998). The small-scale dynamo, as usually 
understood in the literature ( this paper included), is the large- 
Pr m dynamo. 

A smooth random flow is not, strictly speaking, turbulence. 
Turbulent motions at scales above the viscous cutoff are spa- 



tially rough and cannot be approximated by a linear velocity 
field (locally uniform rate of strain). To what extent then is the 
Zel'dovich et al. mechanism applicable to physical systems and 
numerical simulations in which Pr m is not large? 

6.1.1. Suppression of Small-Scale Dynamo at Low Pr m 

A natural way to approach this issue is to consider the oppo- 
site limit: Pr m <C 1 (but Re and Rm are still large). In this case, 
the resistive scale is in the hydrodynamic inertial range and the 
dominant interaction is between velocities and magnetic fields 
at that scale. Is there a dynamo in such a case? 

There currently is no physical understanding of the low-Pr m 
dynamo: the standard heuristic argument involving stretching, 
turbulent diffusion, and resistive diffusion fails to produce an 
unambiguous prediction about the viability of the dynamo be- 
cause all these effects are of the same order and the outcome 
of their competition is a quantitative, rather than a qualitative, 
issue. A number of simulations found in various MHD contexts 
that dynamo shuts down when Pr m is too small (Nordlund et al. 
1992; Brandenburg et al. 1996; Nore et al. 1997; Christensen et 
al. 1999; Maron et al. 2004; Schekochihin et al. 2004a; Haugen 
et al. 2004). In particular, in Schekochihin et al. (2004a), we 
studied this problem using the same code and the same numer- 
ical set up as in this paper and found that dynamo at low Pr m 
is shut down even though the magnetic Reynolds number Rm 
is the same as for successful dynamos with Pr m > 1 . We iden- 
tified two possible interpretations: either a nonhelical small- 
scale dynamo with low Pr m does not exist at all or the criti- 
cal Rm for such a dynamo is much larger than for the case of 
Pr m > 1 and cannot be resolved (cf. Rogachevskii & Kleeorin 
1997; Boldyrev & Cattaneo 2004). 

If a low-Pr m dynamo does exist, what is the structure of the 
magnetic fields generated by it? The Zel'dovich et al. mecha- 
nism does not apply: we cannot have folded fields as there is no 
field-flow scale separation necessary to maintain them. Thus, 
the question is whether the large-Pr m dynamo is the only kind 
of nonhelical small-scale dynamo or there exists yet another, as 
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yet unexplored, variety. 

6.1.2. Small-Scale Dynamo at Pr m ~ 1 

Resolving this issue is important for interpreting the simu- 
lations with Pr m ~ 1. It is clear that the kinematic dynamo 
in such runs is of the large-Pr m kind: we have seen that all 
the main features of the viscosity-dominated runs carry over 
to the case of Re <~ 10 2 and Pr m > 1 (the scale separations 
such as k rms /k\\ change smoothly as a transition from viscosity- 
dominated large-Pr m to turbulent Pr m <~ 1 regime is made: see 
Schekochihin et al. 2004a). This dynamo is driven by the 
viscous-scale eddies and whatever the inertial-range velocities 
might do is guaranteed to happen at a slower rate. However, 
when the stretching component of the viscous-scale eddies is 
nonlinearly suppressed, it is crucial to know if the inertial-range 
velocities that take over can also generate magnetic field at their 
own scale. The heuristic saturation scenarios proposed in § 5.2 
were based on treating these velocities simply as larger eddies 
acting on the magnetic field the same way as the viscous eddies 
did. This description is incomplete if there exists a second kind 
of small-scale dynamo generating magnetic fields in the inertial 
range. In the latter case, the Pr m ~ 1 simulations are in a mixed 
regime containing both the folded fields resulting from large- 
Prm -type dynamo and some other kind of magnetic fluctuations 
generated by the hypothetical inertial-range dynamo (plus, pos- 
sibly, the Alfven-like waves discussed in § 5.2.3). 

The two alternative scenarios of saturation described in § 5.2 
and interpretations of the numerical results in favor of either 
are, thus, qualified by the above discussion. 

6.2. Large-Pv m Dynamo in 3D vs. the 2D Case 

Small-scale magnetic fields in the large-Pr m regime were pre- 
viously studied numerically in two dimensions by Kinney et al. 
(2000). Since the frozen-in property of the field lines holds 
equally well in two and three dimensions, the folded structure 
in thwo dimensions looks very similar to the folded structure 
in three dimensions. Such features as disparity between paral- 
lel and transverse scales of the field variation, power-like tails 
of curvature distributions, and anticorrelation between field 
strength and field-line curvature are present in both cases (Kin- 
ney et al. 2000; Schekochihin et al. 2002b). The reason for 
such similarity is that the folded structure is an ideal effect: 
it is just a geometric property of the field lines arising from 
random stretching by the flow. Diffusion has a selection ef- 
fect on the folded fields with opposite outcomes in three and 
in two dimensions: in three dimensions, the folds can align in 
such a way as to support growth of magnetic energy, — in con- 
trast to two dimensions, where stretching is always overcome 
by diffusion (see § 3.3). Indeed, the magnetic field in two di- 
mensions can be written in terms of a scalar stream function, 
B = (8A/dy,-dA/dx). It is then readily obtained that A satis- 
fies the advection-diffusion equation (Zel'dovich 1957) 

d,A + u- VA = r)AA, (73) 

implying eventual decay of the field. However, the eventuality 
of decay does not preclude transient growth and even a slowly- 
decaying nonlinear "saturation" phase, which was studied in 
some detail by Kinney et al. (2000). They found persistent 
folded structure, near-flat tension spectra, and k~ A spectra of the 
subviscous velocities induced by the tension force — results 
that are analogous to those of §4.1 and §4.2, confirming the 
geometric nature of these features. The key difference is that 
the possibility of sustained dynamo action in three dimensions 



leads to a saturated nonlinear state (§4.3), while in two dimen- 
sions, the nonlinear state is eroded at the resistive timescale and 
the field is destroyed. 

Another feature that arises in three dimensions, but is absent 
in two dimensions, is the possibility for a degree of misalign- 
ment between the reversing fields in the folds: this was mea- 
sured tentatively by A:bj in§3.2.1 and§4.1. We noted in§4.1 
that was larger in the saturated state than during the kine- 
matic growth, perhaps indicating a transition from flux sheets 
to flux ribbons. There is, at present, no quantitative theory of 
three-dimensional magnetic structures and field alignment. We 
have not attempted a full-fledged numerical characterization of 
these features: such a study is left for future work. 

6.3. Diagnosing MHD Turbulence 

In this paper, the set of statistical quantities that were used to 
analyze numerical data was limited to three groups: 

1 . Energy spectra: they show energy distribution over scales 
but miss essential information about the structure of the 
field; also spectra of second-order quantities related to 
the Lorentz force. 

2. Characteristic wavenumbers, pdf's of the field-line cur- 
vature, and correlations between the curvature and the 
field strength — these describe the structure of the field. 

3. The pdf's and moments of the field strength: they char- 
acterize the volume-filling property of the magnetic field. 

All diagnostics in the last two groups (1) are of order higher 
than 2 in magnetic field, (2) contain magnetic field, but not ve- 
locity or its correlations with magnetic field, and (3) are one- 
point in space. That one has to go beyond second-order statis- 
tics is due to the structured and intermittent nature of the mag- 
netic field and to the fact that magnetic back-reaction (B • VB) is 
quadratic in B. Points 2 and 3 are limitations of our description 
of MHD turbulence and deserve some discussion. 

6.3.1. Flo w Sta t is tics 

The results reported in this paper are almost exclusively on 
the statistics of the magnetic field. However, in our discussion 
of back-reaction and saturation in § 4.3 and in § 5.2 we had to 
make conjectures about the way dynamically important mag- 
netic fields modify the flow. The assumptions that helped us to 
construct a saturation model for the viscosity-dominated runs 
were fairly natural: thus, it is clear already from the induc- 
tion equation itself that, in order to stop magnetic field from 
growing, the bb component of the rate-of-strain tensor must 
be suppressed. Indeed, we do find that both (bb : Vu) and 
(|bb : Vu| 2 ) decrease as the transition from the kinematic to 
nonlinear regime occurs. However, in order to resolve the ques- 
tions raised in § 5.2 about the mixing efficiency of the surviving 
motions in the saturated state, the existence of the Alfven waves 
propagating along folds, etc., a detailed and systematic study of 
the velocity statistics and, especially, of the velocity-magnetic 
field correlations, is required. 

Very little numerical evidence on this subject is available in 
the literature. In part, this scarcity is due to the fact that even a 
qualitative physical picture of the magnetic back-reaction has 
been lacking, so it is not obvious what the interesting diag- 
nostics are. Extant results include statements about the growth 
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of (u • B) correlations (Pouquet et al. 1986), positive correla- 
tion between high current and high rate of strain and alignment 
of magnetic field parallel to the vorticity vector (Miller et al. 
1996), suppression of chaos of Lagrangian trajectories by the 
magnetic field (Cattaneo et al. 1996; Zienicke et al. 1998), and 
evidence of nonlocal interaction between large-scale velocities 
and small-scale magnetic fields based on energy-transfer spec- 
tra (Kida et al. 1991; Cho & Vishniac 2000) and on experi- 
ments with switching off velocity field (Maron et al. 2004). All 
these results can be made sense of on the basis of more or less 
straightforward heuristic considerations. While they are all re- 
flections of some underlying back-reaction mechanism, none 
of them tell us explicitly what this mechanism is. Much further 
analysis is necessary to diagnose it correctly. Such an analysis 
is the subject of our ongoing work. 

6.3.2. Two-Point Statistics: Structure Functions 

The second limitation we have imposed on ourselves is that 
all our higher-order statistics are one-point. We believe that 
most of the physics that is qualitatively understood is already 
captured in the diagnostics we employ. Obviously, two-point 
statistics contain a great wealth of information that cannot be 
accessed by the one-point approach. We do not, e.g., attempt to 
study two-point spatial intermittency of the magnetic field. 

In recent years, there has been an explosion of interest in the 
scaling of structure functions in MHD turbulence. These are de- 
fined by S n (y) = (\ip(x+y)-ip(x)\ n ) <~ y^ n , where ip is the quan- 
tity of interest (e.g., the longitudinal velocity field ui = u • y/y). 
The present surge of interest followed the publication by She 
& Leveque (1994) of a phenomenological formula for £„ in hy- 
drodynamic turbulence and the remarkable success it enjoyed in 
reproducing experimental and numerical scalings. Their model 
requires as input parameters the codimension of the most sin- 
gular dissipative structures and the dimensional scaling of the 
cascade (eddy-turnover) time. Applications of their idea to the 
MHD turbulence followed: first, based on the fc" 3 / 2 spectrum 
of Iroshnikov (1964) and Kraichnan (1965) (Grauer et al. 1994; 
Politano & Pouquet 1995), and then on the Kolmogorov k~ 5 ^ 3 
spectrum (Miiller & Biskamp 2000, see Biskamp (2003) for a 
review). The latter model matched numerical data for decaying 
MHD turbulence. In recent simulations, good fits were obtained 
by choosing ad hoc fractional values for the codimension of the 
dissipative structures (Miiller et al. 2003; Haugen et al. 2004). 

Notably, these models seemed to work for structure functions 
not of velocity but of the Elsasser fields z ± = u ± B and it has 
been implicitly or explicitly assumed that magnetic fields and 
velocities are symmetric as in an Alfven wave. However, there 
was no such symmetry in either the decaying isotropic simula- 
tions of Biskamp & Miiller (2000) or the forced ones of Haugen 
et al. (2004). Even in simulations with a mean field, the struc- 
ture functions of u and B, when measured separately, turned out 
to have different scalings, velocity being less intermittent than 
the magnetic field (Cho et al. 2003). On a qualitative level, this 
property has, in fact, been seen in most numerical simulations 
since Meneguzzi et al. (1981). 

Considering Elsasser-field scalings is made attractive by the 
exact result of Politano & Pouquet (1998b) stating that a mixed 
third-order structure function ([z^(x+y)-z T (x)] ■ (y /y)\z ± (x+ 
y)-z ± (x)| 2 ) oc y has a linear inertial-range scaling analogous to 
Kolmogorov's ^ law (this result only holds for Pr m = 1; for ar- 
bitrary Pr m , linear scaling was also proved for two other mixed 
triple correlators — not structure functions — of u and B, see 
Politano & Pouquet 1998a). Thus, while the Politano-Pouquet 



laws provide a set of constraints on u-B (or z + -z~) correlations, 
they do not have exact implications for any structure functions 
involving u, B, z + , or z" alone. 

It is worthwhile to recall that the reason structure functions 
are the diagnostics of choice in hydrodynamic turbulence is the 
Galilean invariance of the velocity field: only velocity differ- 
ences are really interesting. Phenomenological models of inter- 
mittency, the She-Leveque model included, assume that inter- 
actions occur between comparable scales (locally in k space), 
thus giving rise to an energy cascade. Neither Galilean in- 
variance nor A:-space locality of interactions holds for magnetic 
fields: a constant magnetic field cannot be transformed out, and 
an interaction between velocity and magnetic fields at disparate 
scales is exactly what the small-scale dynamo is. The open 
questions are then (1) whether models in the mold of She- 
Leveque are justified for MHD and, if they work, why they do 
so; and (2) whether, in view of the asymmetry between mag- 
netic field and velocity, the Elsasser fields are meaningful. 

Thus, analyzing two-point statistics of MHD turbulence is 
complicated by the unanswered questions about the basic struc- 
ture of this turbulence and by the lack of clarity about what set 
of diagnostics constitutes a physically interesting and informa- 
tive description. We leave these matters to future work. 

7. SUMMARY 

We now itemize the main points of this paper: 

1 . Magnetic fields and turbulence in astrophysical plasmas 
exist and interact at many disparate scales, which makes 
full global description unachievable by brute-force nu- 
merical simulations. We concentrate on small-scale tur- 
bulent fields to probe universal physics independent of 
large-scale object-specific features. In this spirit, we sim- 
ulate incompressible, isotropic, homogeneous, randomly 
forced, nonhelical MHD turbulence in a periodic box. 

2. If no mean field is imposed, all magnetic fields in the 
system are generated by the small-scale dynamo, which 
consists in random stretching and folding of the field 
lines and is best revealed in the limit of large Pr m . In 
its kinematic stage, this dynamo is a one-time-scale pro- 
cess insensitive to the particulars of the random flow. In 
Kolmogorov turbulence, it is controlled by the viscous- 
scale turbulent eddies. It can, therefore, be studied in the 
numerically feasible limit of Re <~ 1 and Pr m > 1 with 
random external forcing. In this viscosity-dominated 
regime, we are able to carry out a parameter scan in Pr m . 
This regime is a model for interaction between magnetic 
fields and the viscous-scale eddies, the forcing represent- 
ing energy input from the larger-scale eddies. 

3. The salient feature of the small-scale dynamo is that it 
produces highly intermittent and structured fields. These 
fields are organized in folds whose length is comparable 
to the flow scale but that contain fields that reverse direc- 
tion at the resistive scale. The scale separation between 
these scales is ~ Pr^ 2 . We diagnose the folded structure 
by measuring various average scale lengths of the field, 
the spectrum of the magnetic tension, the distributions of 
the field-strength and of the field-line curvature, and the 
(anti)correlation between them. The folded structure im- 
plies that the system becomes nonlinear (the fields start 
reacting back on the flow via the Lorentz tension force) 
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when the energy of the folded field becomes comparable 
to the energy of the viscous-scale eddies. 

4. The folded structure is found to persist after the dynamo 
becomes nonlinear and saturates. 

5. Saturation of the dynamo can be achieved via partial two- 
dimensionalization of the velocity-gradient statistics with 
respect to the local direction of the folds, so that satura- 
tion is a result of quasi-2D mixing balancing weakened 
stretching of the field by the anisotropized flow. This 
mechanism is sufficiently robust for a simple semiquanti- 
tative model of saturation based on it to predict magnetic- 
energy spectra in excellent agreement with the numerical 
results for the viscosity-dominated runs. 

6. The folded structure is essential in saturation physics: 
while the magnetic field is formally small-scale (as a re- 
sult of direction reversals), the folds themselves possess 
spatial coherence at the flow scale, which allows them 
to exert a back-reaction on the flow. Both small-scale 
dynamo and this nonlinear back-reaction exemplify the 
nonlocality in k space of the interactions between the ve- 
locity and magnetic fields in MHD turbulence. 

7. While the regime with both large Re and large Pr m is 
not accessible at current resolutions, some progress is 
achieved by simulating the case of Re ^ 1 and Pr m > 1 . 
The small-scale dynamo at Pr m <~ 1 belongs to the same 
class as the large-Pr m dynamo: the magnetic field is gen- 
erated at scales somewhat smaller than viscous and has 
a structure qualitatively similar to the folded structure of 
the fields in large-Pr m viscosity-dominated runs. 

8. In turbulence with large Re, the nature of the saturated 
state depends on how efficiently the anisotropized veloc- 
ity gradients can mix the magnetic-field lines. If mix- 
ing by the inertial-range motions were unsuppressed, the 
magnetic energy would saturate at a below-equipartition 
level comparable to the energy of the viscous-scale ed- 
dies. Numerical evidence, while not definitive, does not 
appear to support such an outcome. 

9. If both stretching and mixing are suppressed we predict 
an intermediate nonlinear regime with magnetic energy 
growing proportionally to time from the viscous-eddy to 
the outer-eddy energy (increase by a factor of Re 1//2 ), 
the fold length increasing to the outer scale (a factor of 
Re 3,/4 ), and the resistive scale increasing via selective de- 
cay by a lesser factor of Re 1,/4 . While these scalings can- 
not be verified at current resolutions, we find clear ev- 
idence of an intermediate nonlinear-growth regime fea- 
turing both fold elongation and selective decay. 

10. The fully developed, forced, isotropic MHD turbulence 
is the saturated state of the small-scale dynamo. While 
phenomenological thinking rooted in the Alfven-wave 
physics of the MHD turbulence with a strong mean field 
leads one to expect a state of scale-by-scale equiparti- 
tion between magnetic and kinetic energy and a Kol- 
mogorov scaling of the spectra, none of the extant nu- 
merical evidence supports this view. Instead, the kinetic- 
energy spectrum is dominated by the outer scale and has 
a steeper-than-Kolmogorov scaling in the inertial range, 



while the magnetic energy is dominated by small scales, 
at which it substantially exceeds kinetic energy. 

1 1 . Based on our numerical results and heuristic arguments, 
we conjecture that the saturated state of the dynamo is 
controlled by the balance of mixing and stretching by the 
anisotropized outer-scale turbulent motions. The mag- 
netic fields in the saturated state are then a superposition 
of folds and Alfven-like waves that can propagate along 
the folds. The existence of these waves in the inertial 
range is predicated on their not being efficient mixers of 
the magnetic fields. These waves are not as yet numeri- 
cally detectable and, therefore, remain hypothetical. 

In the context of our scenarios of saturation with efficient 
and inefficient mixing, it is interesting to recall the two basic 
alternatives put forward more than 50 years ago: magnetic- 
energy saturation at the viscous-eddy energy and scale (Batch- 
elor 1950) or eventual equipartition between the magnetic and 
kinetic energy at all scales (Schliiter & Biermann 1950; Bier- 
mann & Schliiter 1951). Biermann & Schliiter's reasoning as- 
sumed interactions only between velocities and magnetic fields 
at the same scale. Batchelor's argument was based on the anal- 
ogy between the evolution equations for the magnetic field and 
vorticity V x u. While this analogy has since been realized to 
be flawed even for Pr m = 1 (note the recent numerical investiga- 
tions of this issue by Ohkitani 2002; Tsinober & Galanti 2003) 
and the magnetic fields have been shown to propagate to the 
resistive scale, it is still a theoretical and numerical challenge 
to demonstrate that magnetic energy grows above the viscous- 
eddy energy and to explain why it can do so. 

Since the time these alternatives were formulated, theoretical 
predictions and numerical results have ranged from completely 
ruling out the small-scale dynamo (Saffman 1963), to predict- 
ing magnetic-energy pile up at the resistive scale (Vainshtein 
& Cattaneo 1992), to closure-based confirmation of scale-by- 
scale equipartition (Pouquet et al. 1976), to direct numerical 
simulations interpreted as evidence that there is a tendency to- 
ward the equipartition or at least, Kolmogorov scaling (Haugen 
et al. 2003), to direct numerical simulations interpreted as evi- 
dence that such an equipartition is lacking (Maron et al. 2004). 
While our understanding of these alternatives has improved, the 
original dichotomy continues to define the terms of the debate. 

The purpose of the present paper has been to collect a set of 
numerical results and physical considerations emphasizing the 
nonlocal (in k space) nature of the velocity-magnetic field in- 
teraction and the resulting dominance of magnetic fluctuations 
at small scales, which we believe to be the defining properties 
of isotropic MHD turbulence. 
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a In the labels assigned to the runs, "kin" means kinematic regime, "sat" means saturated state. See eq. (60) for the definition of 
Re, eqs. (13-17) for the definitions of the various fc's, eq. (26) for the definition of C k ,b- 

b In the kinematic regime, this run requires higher resolution. 
These runs are slightly underresolved. 
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